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ABSTRACT 


A Finite Element - Boundary Integral Method for Electromagnetic Scattering 

A method that combines the finite element and boundary integral techniques 
for the numerical solution of electromagnetic scattering problems is presented. The 
finite element method is well known for requiring a low order storage and for its 
capability to model inhomogeneous structures. Of particular emphasis in this work 
is the reduction of the storage requirement by terminating the finite element mesh on 
a boundary in a fashion which renders the boundary integrals in convolutional form. 
The fast Fourier transform is then used to evaluate these integrals in a conjugate 
gradient solver, without a need to generate the actual matrix. This method has a 
marked advantage over traditional integral equation approaches with respect to the 
storage requirement of highly inhomogeneous structures. 

Rectangular, circular and ogival mesh termination boundaries are examined for 
two-dimensional scattering. In the case of axially symmetric structures, the bound- 
ary integral matrix storage is reduced by exploiting matrix symmetries and solving 
the resulting system via the conjugate gradient method. In each case, several results 

are presented for various scatterers aimed at validating the method and providing 
an assessment of its capabilities. 

Important in methods incorporating boundary integral equations is the issue of 
internal resonance. A method is implemented for their removal, and is shown to be 
effective in the two-dimensional and Three-dimensional applications. 
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CHAPTER I 


Introduction 


Over the years, integral equation (IE) techniques have been employed for evalu- 
ating the electromagnetic fields scattered from various structures. Special cases of 
the IE methods are surface integral equation (SIE) methods and volume integral 
equation (VIE) methods, which have been the workhorses in computational electro- 
magnetics for the past two decades. They have been used to solve the scattering from 
simple structures, such as smooth perfectly conducting bodies and homogeneous ma- 
terial structures to more complicated ones such as multi-layered coated conductors. 
However, the classical SIE and VIE techniques result in fully populated matrices 
whose required storage is 0(N 2 ), where N is the number of unknowns in the model. 
Consequently, as the size of the structure increases, the memory demand may be- 
come much greater than the storage available on existing computing resources. To 
circumvent this, the IE may be cast in convolutional form by a proper choice of the 
discretization model. The resulting system can then solved via the conjugate gradient 
method, using the fast Fourier transform (FFT) to evaluate the discrete integral op- 
erator. The use of the FFT requires storage on the same order as the dimensionality 
of the scattering body (i.e., 0(Nj ) for three-dimensional scatterers using VIE, where 
Nd denotes the number of cells per dimension), but is still lower than that required 


I 



2 


by the fully populated version. This method, referred to as the CGFFT method 
has been recently employed by several authors to solve two-dimensional (2-D) and 
three-dimensional (3-D) inhomogeneous scattering problems [1, 2, 3]. Nevertheless, 
for 3-D applications the storage requirement remains problematic and furthermore, 
rectangular grid meshing results in a stair-case approximation of otherwise smooth 
scatterers. 

In an effort to further reduce the storage requirement, partial differential equation 
(PDE) methods were considered, since these methods provide for an O(N) storage 
requirement. Among them are the finite difference (FD) methods [4] and finite ele- 
ment (FE) methods [5]. Finite difference techniques typically suffer from the same 
staircasing problems as the CGFFT methods and may have difficulties in modeling 
numerical dispersion when solved in the time domain (TD). Frequency domain so- 
lutions (of primary interest in this study) allow for accurate (conformal) modeling 
of the scattering structure since, for instance, triangles and tetrahedrals can be used 
for discretizing surfaces and volumes, respectively. However, when the finite element 
method (or any other PDE method) is employed, the mesh must be truncated at a 
distance from the structure on which an approximate or exact boundary condition 
is applied. 

Various methods for truncating the FE mesh have been employed. Among them 
are the simple enforcement of an approximate absorbing boundary conditions (ABC) 
at the mesh boundary and the unimoment method. The ABC’s are popular because 
they result in a banded sub-matrix structure. However, they require additional 
unknowns since the enclosure must be placed at a distance approximating the far 
field region. This distance may be large for structures with sharp edges and as a 
result a large number of unknowns is required, especially for 3-D applications. In the 
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unimoment method [6], the scattered field in the unbounded region is represented by 
an eigenfunction expansion. The coefficients of the expansion are then determined 
by enforcing continuity at the circular mesh termination boundary. This method 
produces a dense square sub-matrix whose dimension is proportional to the number 
of modes. It also requires the truncation of an infinite series which may be slowly 
convergent for irregular structures, thus resulting in a large storage requirement. 

An exact termination method is the boundary integral equation, first introduced 
by Hsieh [7] and McDonald and Wexler [8] in the early 1970’s. The method is 
heretofore referred to as the finite element - boundary integral (FE-BI) technique 
and is the subject of this thesis. The boundary integral equation employs the free- 
space Green’s function which implicitly satisfies the radiation condition at infinity. 
Since the integral equation is exact, the mesh termination boundary may be brought 
very close to the scatterer to minimize the meshing of the free-space regions enclosed 
by the termination boundary. The main drawback of the method is again the dense 
submatrix structure associated with the boundary integral equation. However, by 
carefully choosing the shape of the boundary, some of the boundary integral terms 
are cast into convolution form. Thus, on employing a conjugate gradient solver, 
these integral operators may be evaluated via the FFT as was done in the CGFFT 
method. The required storage of the boundary integral is then reduced toward 0(N) 
and allows for the accurate solution of large inhomogeneous scatterers as well as small 
targets. 

In this thesis, the FE-BI technique (presented in general terms in chapter II) is de- 
veloped for two-dimensional and axial-symmetric structures. The two-dimensional 
case is based on [9]. Chapter III contains a FE-BI formulation for rectangular enclo- 
sures, leading to simple boundary integrals some of which have convolutional form, 



4 


thus leading to a reduction in the memory requirement. In chapter IV, circular and 
ogival enclosures are considered. For the circular boundary, the boundary integral is 
entirely convolutional in form and an O(N) storage requirement is thus achieved at 
all times. Circular enclosures are consequently preferred for storage efficiency if the 
structure’s outer boundary does not substantially deviate from a circle. 

In chapter VI, the FE-BI formulation is developed for axially symmetric struc- 
tures. The finite element method for this problem was originally developed by Mor- 
gan and Mei using the coupled azimuthal potential (CAP) equations for generating 
the FE matrix system. In that formulation, the unimoment method was employed for 
terminating the mesh, whereas in our implementation the boundary integral equation 
is instead employed for this purpose. It is related to Flemings implementation [10] 
for circular boundaries, but the present formulation allows for an arbitrarily shaped 
enclosure. Consequently the boundary integrals are no longer convolutional and the 
storage reduction is in this case achieved by exploiting certain symmetry properties 
of the boundary integral subsystem. A problem with the CAP equations is the pres- 
ence of a singularity in the finite elements which tends to corrupt the solution, and 
a method is presented for circumventing this difficulty. 

A difficulty with most boundary integral formulations is the appearance of in- 
ternal resonances which corrupt the solution. These resonances correspond to the 
eigenvalues of the integral operator, and also occur at the same location as the cut- 
off frequencies of a resonator with conducting walls of the same shape as the mesh 
termination boundary and filled with the same material as the unbounded medium. 
The resonances are particularly problematic when the structure becomes electrically 
large, in which case the eigenvalues are very closely spaced. Without any correction, 
computations for large structures are unreliable, and a method for correcting the 
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resonance problem is developed in chapter V. 

Specific contributions of this work include the implementation of the FE-BI for 
rectangular, circular and ogival boundaries as described in chapters III and IV in 
a manner taking advantage of storage reduction schemes. A new method is also 
presented for suppressing the resonance corruption problem existing in almost all 
implementations employing some form of a boundary integral equation over a closed 
surface or contour. This method involves the introduction of a complex wavenumber 
and is demonstrated for both 2-D and axial-symmetric bodies. The implementation 
of the FE-BI method for axially symmetric structures as described in chapter VI 
is for the most part new and incorporates a scheme for treating the line singularity 
associated with the CAP equations. This resulted in real matrices for lossless scat- 
tered, a property consistent with 2-D and 3-D implementations of the finite element 
method. 



CHAPTER II 


Fundamental Concepts 


In this chapter, some basic concepts are presented as applied to electromagnetic 
scattering and its computation via the FE-BI method. At the end of the chapter, 
the conjugate gradient matrix solver is also appended. 

2.1 Basic Electromagnetic Theory 

The overall goal in any electromagnetic problem is the solution of Maxwell’s 
equations subject to given boundary conditions. In the frequency domain, Maxwell’s 


equations in free space are 

Vx£= —jufiH (2.1) 

V x H = jutE (2.2) 

V • (eE) = 0 (2.3) 

V • {(iH) = 0 (2.4) 

where 

E(r) — electric field strength [Volts/meter] (2.5) 

H(r) = magnetic field strength [Amperes/meter] (2.6) 

e(r) = electric permittivity [Henrys/meter] (2.7) 
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p(r) = magnetic permeability [Farads/meter] (2.8) 

Note that the last two equations, which are based on Gauss’ law, can be obtained by 
taking the divergence of the first two equations. Thus for time harmonic fields, only 
the first two of Maxwell’s equations are needed for a unique solution of E and H . 
By substituting (2.2) into (2.1) and vice versa, we obtain 

V x —V x E — u> 2 tE = 0 (2.9) 

I* 

V x jV x H - uj 2 pH - 0 (2.10) 

which are commonly referred to as the vector wave equations. They must be solved 
and are subject to the boundary conditions on the particular scatterer as well as 
the radiation condition. Typically when solving (2.9) and (2.10) for a piecewise 
homogeneous scatterer, the following conditions must be explicitly imposed: 

• Continuity of tangential E and H 

• Continuity of normal pH and tE 

• Tangential E = 0 on metallic surfaces 

For radiation and scattering problems, the radiation condition must be also satisfied. 
This condition assumes that the wave must be outwardly propagating and must 
attenuate no slower than the inverse of the distance far from the source. Expressed 
in mathematical terms, 

jinn [tjH — r x E] = 0 (2.11) 

him [E - f x rjH] = 0 (2.12) 

where tj = is the impedance of unbounded medium and f is the unit radial vector 


in spherical coordinates. 
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The media considered is this work is linear and isotropic and the material pa- 
rameters fx and e are often normalized to their free space (vacuum) values of e 0 — 
8.85 x 1CT 12 H/m and fi 0 = 4 tt x 10" 7 F/m. The relative constitutive parameters of 
the medium will be denoted by 

e r = - (2.13) 

fr = f (2.14) 

Mo 

A fundamental solution of the wave equation is the plane wave 

E(r) = t~ fkT (2.15) 

where k = Uy/jlt and u = 2irf is the frequency in rad/sec. For lossy media, // 
and e are complex and thus k = k T + j'fc,-, implying an exponentially decaying field. 
Note that the required condition, ki < 0, ensures that the wave does not grow 
exponentially and, consequently, the imaginary part of n and t must likewise be less 
than zero. 

2.2 The FE-BI Approach 


To solve PDE’s (2.9) and (2.10), the scattering body is first enclosed in an artificial 
surface S bounding a volume V and for this application all sources will be assumed 
to reside in the unbounded region outside S. To satisfy Maxwell’s equations in this 
volume, the vector wave equation for the electric field in (2.9) is discretized via the 
method of weighted residuals. Though not given here, a similar approach in terms 
of the magnetic field may developed in a parallel fashion. 

The volume V is first divided into N e elements. Within each element, the 
weighted residual expression is given by 



- V x E - u?tE)dV = 0 
M 


(2.16) 
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where W‘ is the ith weighting function over the eth element. First, we temporarily 
define A as 


and recall the identity 


A = 


-V x E 


(2.17) 


W' • V x A = V • (A x W { ) + A • V x W- 


(2.18) 


Using this identity, the first term in (2.16) can be rewritten as 



• V x AdV e 


JJJv-(Ax W])dV e + /// 4 • V x 1 V'dX 
V e V e 

fjn e -(Ax W])dS e + JJJa-Vx W]dV e 

S' 


(2.19) 


where the last equality was obtained by invoking the divergence theorem and thus, 
S e is the surface enclosing V e and n e is its outward unit normal. Substituting (2.19) 
into (2.16) yields 


III (-Vx£)-(Vx W\)-u?tE 

ye L V 

and upon rearranging this we have 


dV e 




uPtE 


dV e 


juW] ■ (n e x H)dS e 
S e 


( 2 . 21 ) 


For scatterers comprised in part of conducting material, the surface boundary of 
the volume V is S = S a + S c , where the subscript a denotes the exterior boundary 
and the subscript c denotes the conducting boundary. In the assembly of the finite 
element equations, the fully discretized form of (2.21) is summed over all elements. 
Since the tangential fields are continuous across material interfaces, the surface inte- 
gral in (2.21) will cancel everywhere except on S a and S c . Thus, unless S e intersects 
S a or 5 C , the surface integral may be removed. 
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To construct a system of equations from (2.21), the field in element e may be 
represented as 

I? = f^W j (x,y,z)E e j ( 2 . 22 ) 

j = i 

m 

n x ZT = xl\F j (x,y,z)E* e (2.23) 

3=1 

m 

nxF = ^nx 7V'(x,y, 2 )//f (2.24) 

j=x 

where the superscript t denotes the tangent to the boundary of element e. (These 
vector shape functions must also be chosen to satisfy the divergence condition given 
by (2.3).) Substituting this expansion into (2.21) and using Galerkin’s approach (i.e., 
W* = TV*) gives 

[< v x *1) • < v x dV ‘ 

3 = 1 y e " 

= E ff ■ (« # X + E $ 3»N*i ■ X A^)d5 e (2.25) 
j=1 5 e € 5 a J =1 5'eS c 

where A: = Uy/JIi. 

After summing over all elements, our system takes the form 

{e:} 

{ed 

{Ei} 

{E?} 

{ED 

m 

In these the subscript a refers to unknowns on iS a , c refers to the set of unknowns on 
S c and / refers to the remaining unknowns. Furthermore, the superscripts n and t 
refer to the norma] and tangential components, respectively, of the weighting and/or 



at; k, k ;< at; o 
A'“ A-“ K l al K" /C s. 

AT. A'i. A„ Kj c I<; c 0 

AT." A"' A' t ", A" n A”' 0 

A'“ AT AT, A'T AT 0 
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testing functions. The application of the boundary condition on S c results in the 
following system 

{£„"} 

{£■;) 

{£/} 

{£,") 

{£;} 

{HD 

Clearly, another equation is necessary relating E l a and H[ on S a - The FE-BI 
approach employs the Stratton-Chu integral equation for the purpose. This equation 
relates the electric and magnetic fields via the integral equation 

E(r) = T(r) + ff {-j* 0 [n' x r 1o H(f , )]g(r, 7) + jj-[n' • V' x T /o 77(f')] V'g(r, ?) 

+ [n' x £( 7 ^)] x 

(2.28) 

where r/ 0 = \J fi 0 A 0 and g is the free space Green’s function. Evaluating this equation 
on the surface S a and discretizing it, we obtain the subsystem 

L' aa {K} + M< aa {K} = U' a (2.29) 


Ar k% i<:, k% o o 

K: Ka Ki A :« 0 B a 

K?a K\a A I<\ c 0 0 
A'r K2 IQ K™ 0 0 

0 0 0 0 /0 



which relates the tangential fields on S a . Combining (2.29) and (2.27) we obtain the 
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block matrix 


K nn 

aa 

L'nt 
A aa 

K, 

Kl 

Rtn 

aa 

K aa 

Kl 

K“ 

K n 

n Ia 

Ka 

Kn 

K\ c 

I<™ 

ca 

A '" 1 

ca 

K 

Kc 

0 

0 

0 

0 

0 

LL 

0 

0 


0 0 

0 B aa 
0 0 
0 0 

1 0 

0 ML 



{£.") 


0 


(£■;} 


0 


{£/) 


0 


{*?} 


0 


(£;) 


0 


m 


. {u ‘ ] . 


(2.30) 


which, after applying the condition on 


S c , can be solved uniquely for the interior 


and boundary fields. The resulting system is then solved via the conjugate gradient 
method. 


After the fields in the solution region have been determined, the scattered fields 
are computed by evaluating the integral equation in the far zone. In particular a 
quantity called the echowidth is computed for 2-D problems and is given by [11] 


a = lim 2 np 

p—*oo 


l^’l 2 


(2.31) 


where p is the usual cylindrical coordinate and <f> denotes either E z or H z . For three- 
dimensional applications, the quantity of interest is the radar cross section (RCS) 
and its expression is 


a = lim 47rr 2 


|£-(r)p 


(2.32) 


where r is the radial spherical coordinate. Expressed in component form we have 


a pq = lim Airr 2 — - 

r— ►co 


|£‘(r)| 2 


(2.33) 


where p and q represent the polarization of the scattered and incident fields, respec- 
tively. 
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2.3 The Conjugate Gradient Method 

The CG algorithm is employed throughout this work and is given as follows for 
solving the system Ax = b. 

Initialize the residual and search vectors 

76 = II [ 4 > inc o o o ] T IlHl b Hi 

5 = A<j> {0) 
r* 1 * = b — s 
s = A a r [l) 

7, = II - \\l 
/? ( °) = 7 ;i 

p W = 0(O) S 

Iterate for A: = 1, N g 

s = Ap w 
Is = \\s \\l 
« (fe) = 77 1 

^ k+ 1 ) = <£(*) + Q ( k )pW 
r (*+i) _ r (k) _ Q (k) p (k) 

lr = II \\l 

S = A a r ( * +1) 

Is = || - \\l 
p (k) = 77 1 


P ( fc+ x > = pW + 
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Terminate when k = N g or < tolerance. 

The given tolerance must be used with caution. In each of the FE-BI techniques 
developed, the excitation b is only partially full. If the simple preconditioner of divid- 
ing by the diagonal value is employed, then those rows associated with an excitation 
component of zero will consequently not scale the excitation. In other words for 
those rows, 6 will not change in value while the matrix itself does. As a result, the 
residual error will be scaled differently and, consequently, the tolerance criterion will 
be scaled differently as well. This may lead to artificially small normalized resid- 
ual errors and cause the convergence criterion to be satisfied before convergence is 
actually achieved. 

One way to avoid this is to compute the far field in terms of a at one or more 
angles every few iterations and compare the current value to the previous one. When 
the difference is sufficiently small for 15 or 20 consecutive iterations, convergence is 
assumed to have been achieved. 



CHAPTER III 


A Finite Element — Boundary Integral 
Formulation for Two-dimensional Scattering with 
a Rectangular Termination Boundary 

3.1 Introduction 

In this section an FE-BI method is developed for two-dimensional scattering in 
which the exterior boundary is rectangular in shape. This enclosure ensures some of 
the boundary integral terms are convolutional and are therefore amenable to evalu- 
ation via the FFT in the conjugate gradient solver. Results are presented for several 
structures. 

3.2 Analysis 

Consider a cylindrical body of arbitrary cross-section and composition illuminated 
by the plane wave 

J nc {p) = z^ip) = ze’ koPCOs{e - 6 ^ (3.1) 

as indicated in Fig. 3.1. To evaluate the fields scattered from this object, two 
boundaries are placed tightly around the body as shown in Fig. 3.2. Inside the outer 
boundary, the Finite Element Method is applied to solve the Helmholtz equation 
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Figure 3.1: Geometry of the scatterer 



H Ah 


Figure 3.2: Partially discretized body 


H> f- 


17 


given by 


V • [v{p)V <t>{p)] + k 2 0 u{p)<f>{p) = 0 


where 


for E-polarization and 


<Kp) 

v{p) 

u{p) 


E Z {P) 

1 

Pr(p) 

trip) 


m = h 2 (p) 

v(p) = —TTT 

trip) 

“(p) = Pr{p) 


(3.2) 

(3-3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 


for H-polarization. Also, k 0 = ujy/p^tl is the wave number, and p r and e r are the 
relative permeablility and permittivity, respectively. 

The appropriate boundary condition is enforced on the surface of the impenetra- 
ble boundary, while the radiation condition is satisfied implicitly by evaluating the 
integral equation 


m = r c (p) - / rt |g(?, 

on the boundary T a , where 


p) 




-m 




V 


(3.9) 


(3.10) 

is the 2-D Green’s function in which H^\') denotes the zeroth order Hankel function 
of the second kind. Furthermore, ~p and ~ff are the usual observation and source 
position vectors, respectively, and 


|p-p| = - x ') 2 + {y - y ') 2 


(3.11) 
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Symbol 

Description 

N g 

number of nodes in the finite element mesh 

N e 

number of elements in the finite element mesh 

N x 

number of nodes on T a or Tt> along the avdirection 

N y 

number of nodes on T 0 or Tt along the y-direction 

N a 

total number of nodes on T a 

N b 

total number of nodes on Tj, 

N ab 

Na + N b 

r a 

ELi r 0j 

V 

ELi n. 

r c 

E?=i r c , 

(^a,} J/a,) 

coordinates of a point on contour r ai 

(*6,» Vbi) 

coordinates of a point on contour Tj,, 


Table 3.1: Definitions for the finite element mesh 
The normal derivatives are taken in the direction of the outward normal of r c . 

3.2.1 Discretization of the Object and Field Quantities 

In Fig. 3.2, T a is the field/observation point boundary, and T c is the integration 
contour, which is placed midway between T a and IV Also, Tj denotes the contour 
enclosing the impenetrable portion of the scatterer. Herewith, each side of T a , r 6 
and r c are numbered counterclockwise starting from the top side, as indicated in 
Fig. 3.2. The fields in the region between T a and Tj satisfy (3.2) in conjunction with 
the required boundary condition on Tj. The boundary integral equation (3.9) will 
be enforced on T a . 
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Symbol 

Description 

K 

nodal fields on T ai 

<t> b, 

nodal fields on Tt, 

K 

nodal fields on 

4>d 

nodal fields on Tj 

4>i 

region enclosed by r& and Tj, exclusive 


Table 3.2: Definitions of the field vectors 

To numerically solve (3.2), it is required that the region within T a be discretized. 
This is done in a traditional manner when employing the finite element method. 
The global node numbering starts from the right endpoint of contour T ai and pro- 
ceeds counterclockwise. The numbering continues beginning at the right endpoint of 
contour T^ and proceeds counterclockwise. Within IT, the nodes are numbered arbi- 
trarily. The definitions pertaining to the FE mesh are given in Table 3.1. Each node 
corresponds to an unknown field value to be determined. It is important to associate 
the unknown field values corresponding to the various node groups on contours T a 
and IT by using different variables. The labeling scheme is given in Table 3.2, and 
this discrimination of the nodal fields is required, since they are treated differently 
in the analysis. 

3.2.2 Derivation of the Finite Element Matrix 

One of several approaches may be used to generate the finite element matrix, such 
as the variational approach or the method of weighted residuals. In this development, 
we will utilize Galerkin s method, which is a special case of the method of weighted 
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residuals with the distinction that the testing functions are the same as the basis 
functions. 

Proceeding with the finite element analysis, we may rewrite (3.2) as 

^ ^ ^ ’ Q 

!h + ^ v ( x ,y)fojH x ,y) + *%u(x,y)4{x,y) = 0 (3.12) 

the residual of which is given by 

d m 0 ’ 0 ’ 0 " 

R = ~]h v ( x 'V)fctt x *v') ~ v ( x ^y)^( x ^y) - k lui x ,y)<K x ,y) (3.13) 

The field within T a may be represented as a summation of piecewise continuous 
functions and, thus, may be written as 

<t>{ x ,y) = ^2<f> e i x ,y) (3.14) 

c=l 

where <j> e (x,y) is the field within element e. It is expanded as 

<i> e {x,y) ^ it N j( x ’y) $ (3-15) 

3 = 1 

where Nf’s are the standard shape functions (found in any standard finite element 
book), tf's are the fields at the nodes, and n is the number of nodes per element 
(n = 3 for linear elements considered here). The weighted residual equation for the 
eth element is defined by 
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Further, by using the identities 
d ( dN'\ 

N 'r x hr 


d 


dx V dx 


dN ' 

J N' 


— v 


dNf dN] 

dx dx 




dN. 


dN' dN' 


-u- ‘ J 


dy \ dy ' J dy dy 

and the divergence theorem 

// (t; + %) dxdy = L ft ' {ui + vy) dl 

where C' is the contour enclosing the eth element, (3.17) becomes 

" ff ( dN' dN' dN' dN' \ 

£ // + v ~dy'~at " k ° uN ' N ’) dxdy 

n r ( dN' dN' \ 

= ^/ ct + v ~dy~^J dl 1 = 1 ’ 2 ’ -’ n 


:=l S' 


This may be written in matrix form as 


(3.18) 

(3.19) 


(3.20) 


(3.21) 


A'(f>' = b' 


where 


n- II (•%%>•%■■ ?-«»:)« 


and 


n r ( dN' dN' \ 

{ fc '}. = E£ e ^ e + -” rf/ * = 1 > 2 — ’ n 


For linear triangular elements, N' are given by 


1 


N' = ^7 K + b'x + c'y) 


and 


Cl' = - 


i y\ 

1 x\ y\ 
1 x% y' 2 


= M } ~ b'c') 


(3.22) 


(3.23) 


(3.24) 


(3.25) 


(3.26) 
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< = x )vl - xty- 

(3.27) 

S* 

1 

II 

v ■« 

-o 

(3.28) 

c' = x e k -x) 

(3.29) 

where (x-,y, e ) is the coordinate of the ith node of the eth element. 

From (3.25) 

dNf __ b' 

dx 2U e 

dN* c* 

dy ~ 2f V 

(3.30) 

(3.31) 

Substituting these and the formula [5] 


JjWYWY ixty = 2 

§* (P + <7 + 2)! 

(3.32) 

into (3.23), we find 


\ A % = 4Q. m + «<$) - i + ««> 

(3.33) 

where 


*,={ iifi=> 

1 0 otherwise 

(3.34) 


In (3.33) we have assumed that u and v (the material constitutive parameters) are 
constant within each element and are given by u e and v e , respectively. By summing 
over all elements as implied by (3.14), we may write the overall system in block form 


A aa A a b 0 0 


<f>a 


b a 

Aba Abb Abl 0 


<t>b 


0 

0 Aib An Aid 


<t>I 


0 

0 0 Adi Add 


J>d _ 


0 


(3.35) 
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In this, the values of the elements in the submatrix A vq are the contributions asso- 
ciated with the nodes in group p which are connected directly to the nodes in group 
9- 

One can easily show that the line integral contribution (3.2-1) of those elements 
vanishes everywhere, unless the element has a side in common with T a . As a result, 
6 e contributes only from the boundary r a of the finite element region, as indicated 
by the presence of the vector b a in (3.35). Without a priori knowledge of the total 
field on that boundary, 6 “ cannot be determined. We may, however, provide the 
appropriate condition on this boundary by utilizing the integral equation (3.9). This 
amounts to replacing the first block-row of the matrix (associated with 4> a on T 0 ) 
with a discrete form of this integral equation. 

3.2.3 Evaluation of Boundary Integral 

The boundary integral in (3.9) may be written as a summation of four integrals, 
one for each side of the contour T c as 

■/r„ Gtoft-sr*#) + dl n 

~l„ di a 

- j r ) f \ dl c , (3.36) 

where the derivatives along the x and y directions, denoted by and ab re ‘ 
spectively, have been left in this form for the later convenience of determining them 
numerically. More explicitly, we have 


<t>(x,y ai ) = 4> inc {x,y ai ) 

3 3 
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“ J r |g(x - x', y 0 j , y Cl )-^—<j)(x', y Cj ) — <f>(x\y Cl )-^G(x - x\y a ^y Cl ) 

f \ 9 9 

G(x, x C2 , y a , ,y <f>{x C2 ,y ) 4" ^(xejiJ/ ) ^*(x, x C2 ? 2 /a i •* J/ ) 


dx' 


dy' 


I G(x x,y ai , y C3 ) <^(x , y C3 ) ^(x , y C3 ) ^ G(x x , j/ a j , y C3 ) 

- r 


afx' 


<9 <9 

G(l, X Ci , T/o i i J/ ) ^(*C4 > 2/ ) 4*{ x ci j 2/ ) t^^G(x, x C4 , J/o j , y ) 


dy' 


(3.37) 


and 


<f>i x a 3 ,y) = 4>' nC ( x a J ,y) 


9 9 

G(x a 2 , x , j/, t/ Cl , y Cl ) — , y C j )^^ 7^( :Ca 2 ? ^ * !/? J/ci ) 


G(x a j , ^C 2 i y y y ) "t~ ^(^C 2 ? y ’ *^ c 2 ’ y y ) 


-I, 

- 1, 

^ 1 1?! X a 2 ? ^ } JA J/ca)^^ * J/ca) 4" i ^ c 3 ) 2 1 % yV> Vc^) 

-l 


dx' 


dy' 


dx' 


9 9 

G(x a ,,x Ci ,y - y')-^—<}>{x c< ,y') - <f>(x Ct ,y')—G{x a ^x Ci ,y- y') 


dy' 

(3.38) 


where the first subscript on x or y refers to the contour (a, b or c), and the second 
refers to the contour number (see Fig. 3.2 and Table 3.1). It is noted that the 
arguments of the Green’s functions have been modified to imply a convolution when 
appropriate, and this representation will be used throughout. 

The normal derivatives of <f> may be evaluated via the central difference formulas 

9 1 

C y ') = ^ Mx„, y') - <j>{x h , y')] + 0( A 2 ) (3.39) 

9 1 

Q— <K X \ Vc) = ^ [4>(x', y a ) - y fc )] + 0( A 2 ) (3.40) 

where A is the displacement of T 0 from r 6 (A is usually less than one tenth of a 
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wavelength). Substituting (3.39) and (3.40) into (3.37) and (3.38), we obtain 


) = <(> tnc (x,y ai ) 

3 3 

- j \l< ~ G{x - x',y ai ,y Cl )(f>(x’,y ai ) - K+G{x - x',y a , ,y Cl )(f>(x',y bl ) 

*/r Cl L 3 3 

- I \K^G{x,x C2 ,y a] ,y')<f>{x a2 ,y') - K~G{x,x C2 ,y ai ,y')<f>{x b2 ,y') dy' 

Jr C2 l 3 3 J 

(7(x x , y a , , y C3 )(f>{x , y a3 ) I\y G(x x , y a , , y C3 )d>{x , t/i> 3 ) 

- I \K~G(x,x Ci ,y ai ,y')<t>(x at ,y') - A'+G(x,x C4 ,y a , ,y')4>{x b ^y') dy' 

Jr CA L 3 3 J 


dx' 


dx' 


(3.41) 


and 


i y) $ (*^02 9 y ) 

4 4 

- [/f“G(x 0 , , x', y, y Cl )</>(x', y 0 j ) - /f J G(x a 2 , x', y, y Cl )<?!>(x', y 6l )] dx' 

- L h +G ‘ x a2 ,x C2 ,y - y')(f>{x a2 ,y') - K x G{x a2 ,x C2 ,y - y')4>(x b2 , y')j dy' 


dx 


- L h +G < 2 1 ^ 1 y 'I J/c3 )^(*£ 1 J/03 ) Ky G(x a2 ,X ,y,y C3 )<^(x 1 J/63 ) J 
/ l^x 2 ? ^ y y ) 4 > {^'ai^y ) ~~ K x G(x a 2 1 ^C4 1 y — y ) 4 > {^bi^y ) dy 

d r c< L 4 4 

(3.42) 


in which 


A -± _ J_ , IA 

* A 23x' 

,.± = I , I A 

» A 2 dy 1 


(3.43) 

(3.44) 


Assuming a pulse basis expansion for the nodal fields (i.e., piecewise constant func- 
tions centered at nodes on contours r o and Tj,), a midpoint integration may be 
performed for the evaluation of the integrals in (3.41) and (3.42), to obtain 


<t>{xi,Va , ) = <f> ,nc (xi,y ai ) 
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Nr r 


r'y ^j^Vai ^yci) 4 > { x jiVaj) Ry G(x, Xj , y a j , y Cl )^>(Xj , l/jjj ) 

J=1 L 3 3 

r 

~ ^ XI A 'x @( x ii X C2 > J/o j i »j)^(*« a i J/j) A r G(x,, x C2 , y a , ,y J )(f>[xf >2 , ?/_, ) 


j=i 

W. r 


- A XZ [ A 'y + G (*. - Xj , y Q , , y C3 )<f>(xj,y a3 ) - A' y G(x, - x i? y a , ' J/c 3 )<^(x J , J/(, 3 ) 

% r 

~ ^ XZ A x G(x{, x Ci ,y a , , yj)<f>(x at ,yj) — K+G(xi, x Ct , y a , , yj)<f>(x-b t , t/j ) 

j=l L 3 3 , 


(3.45) 


and 


<£(*a 2 ,!/,) = <f> ,nC (x a2 , J/,) 

4 4 

N* r 

~~ ^ XI Ay ^(' r °2 ? X j> Vi) »c, yaj ) A* ^*(^02 , Xj, ?/i, J/cj )^>(Xj, ) 


3 = 1 L 


- AEk + G(x« 3 , 

i-i L 


JV« 


A XZ 

i=i 


/CG(x„,. 



X c 2 , J/t yi)^(^a 2 ^ J/j) A x G(x a2 , X C2 , y t — J/ ; )0(xt 2 , J/j) 

4 

yiiVcs ) ( t > ( x j * J/a 3 ) — A y G(x a2 , Xj, y,, y C3 )<^(xj, yt 3 ) 
x c 4 , J/i “ J/i)^a 4 ,I/i) - A^G , (x a2 ,x C 4 ,y t - t/j)<£(x 6 4 ,yj) 

4 


(3.46) 


In these x,- and y,- denote the ith matching/ testing points corresponding to the nodal 
locations on r a , while Xj and j/j denote locations on Tj,. We recognize some of the 
terms in (3.45) and (3.46) as discrete convolutions amenable to numerical evaluation 
via FFT. The subsystem (3.45) and (3.46) may be written more concisely as 


<t>a 1 


6 inc 

* a \ 


6 , 

7 1 + 

s t, h R* 

^ 0 , 6 4 


^a, 

4*a 2 


6 inc 
r a 2 


T a 2 b l 



■S’ajfc* A y 


<&» 2 

<Pa 3 


6 inc 

’ 03 


s: 3bl Rr 

rife 

c+ 

^<*363 

^64 


^<>3 

0«4 _ 


6 inc 
' a 4 


_ T Zb , 


T+ 

a 4 fc 3 

5 ’a 4 6 4 


^a 4 
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si* 

T~ 

a ih 


K* 

-o 
-U ^ 

^ 'a 2 6 2 

T~ 

a 2 6 3 

Si*R, 

si*R, 

T~ 

J a 3 6 2 

^ 3^3 

K* 

Km 

^462 Ry 

0463 

Km 


(3.47) 


with the various parameters to be given explicitly later. The matrices R iy simply 
reverse the order of the unknown vector so that the convolutions may be performed 
properly. This is required solely because of the employed counterclockwise nodal 
numbering scheme. 

Since 


the vector 



i = 1,2, 3, 4 

j = 2,3,4, 1 


-.T 


^6, 0t> 2 ^63 ^fc 4 


(3.48) 


(3.49) 


can be related to the actual unknowns on the contour Tj, via a transformation D b as 


4>'b — Db<t>b 


and (3.47) may then be written as 


(7 + L aa )4> a — L a bDb<t>b = 4>'a C 


(3.50) 


(3.51) 


or 



<Pb 


j tnc 
= <Pa 


(3.52) 
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where 


/ + L aa — 


Lab — 


I + S a i b 1 

T+ 

I a 1 b 2 

C + 

°ai6 3 

R x 



T~ 

aib i 

I + SU T+ b3 


S a 2 b i J 

Ry 

^a 3 6i 

T+ 

Q362 

1 + si.. 

T~ 


T„t, 

S* lb ,R, 

T .\ 


I + S 

a 4 6 4 


T~ 

<*ih 

F ai b 3 F X 

• K 

b* 


T+ 

<*261 

^a 2 b 2 

T~ 

02&3 

st 

64 Liy 


^0361 R* 

-^0365 

S a 3 b 3 

K 

64 


K<n 

Sa,b 2 F y 

T a<b 3 

st 

bi 



(3.53) 


(3.54) 

After replacing the first block row of (3.35) with (3.51), the complete system may be 
finally written as 


/ + L aa 

LabD b 

0 

0 

Aba 

Mb 

Mi 

0 

0 

A/6 

A n 

Aid 

0 

0 

Adi 

Add 


(3.55) 


to be solved via the CG algorithm. 

The elements of A BI defined above may be evaluated via the discrete Fourier 
transform. Specifically, we have 

= DFT ~ 1 { DFT [G(x,y ai ,y bl )±G y (x,y ai ,y bi )}DFT[cf> h ,]} (3.56) 

S ti > >> = DFT ~ l { DF T[G(x,y a3 ,y bl J±G y (x,y a3 ,y bl ))DFT[<t> bl J ] i (3.57) 

= DFT-' {DFT[G(x aj ,x b] ,y)±G x (x aj ,x br y)}DFT[<j> b ,} ] j (3.58) 
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= DFT 1 ^DFT[G(x ai ,x h z,y) ± G x {x ai ,x b2 ,y)]DFT{<t> b2 }^ (3.59) 

in which DFT denotes the discrete Fourier transform operator. Also 

G(x,x',y,y') = -j-H {2) (k 0 \p - p'\) = -^H {2] (k oy f[x. - x') 2 + (y - y') 2 ) (3-60) 

G x (x,x',y,y') = ^-^pG(x,x\y,y') 

j . . H[ 3) (k 0 y/(x - x') 2 + (y- S/') 2 ) . 

= -£AA; 0 \ V M* - x ) (3.61) 

8 y/{x - x') 2 + (y - V ) 2 

Gy{x,x',y,y') = y^ 7 G(x,x , ,y,^/ , ) 

,3.62) 

8 \J (x - x ') 2 + {y- j /') 2 

Special cases of the convolution operators for the chosen mesh are given as 


G x (x a2 ,£/5 2 , y y ) — 


j ni 2> (hj(77, -n, + <!/ - sot) 

® " - 14, ) 2 + (!/ - j') 2 


- x t - |A 


( 3 . 63 ) 


Gy (x x , $/„ j , j ) — 


*8*’ 


/(x - x') ! + (y„, - y tl y 


■|y«i - yt. I A 


( 3 . 64 ) 


and the corresponding expressions for G are implied by the arguments in the previous 
two equations. Additionally, the upper coordinate in the braces corresponds to the 
upper sign and likewise for the lower one. 

The cross-term element submatrices are given by 


Ta 162 G^Xi) iVai ? Vj ) i G *^62 1 2/ a 1 1 Uj ) 

3 4-1 ij 4 3 4 3 

2 6j 2 1 1 J/t? J/6 j ) i ^y(*^a2 ** ^ j ? J/n 2/6 1 ) 

4 3 J 1 1 4 3 4 3 


( 3 . 65 ) 

( 3 . 66 ) 
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G x {xi, Xb 7 , y a j , j/j ) — 


• ( k ° y(x,-x fc J 2 + ( y a -Jo) 2 ) 

=Fj*o 7 V ki - x b2 |A 

8 J(xi - ar >2 ) 2 + {y ai - y ,-) 2 4 


(3.67) 


Gy (X a 2 5 3/»1 3/6 J ) 


r ■ - * ■ 3 \ yi - y bi |A 

8 ,/( X “! “ X j ) 2 + ( Vi ~ Vby ) 2 3 


(3.68) 


where again the corresponding expressions for G are implied by the arguments of the 
previous two equations. Making the substitutions 


(x,-*o 2 = (;-!) 2 a 2 


(3.69) 


(y ai — Vj ) 2 = j 2 A 2 

/(x» - Xfcj ) 2 + (y a , - yj) 2 = A\/(* - |) 2 


(3.70) 

(3.71) 


(x 02 -x > ) 2 =i 2 A 2 

4 

(y. - y», ) 2 = (» - |) 2 a 2 

3 

/(x a2 - Xj ) 2 + (y, - y hl y = At // 2 + (i - 


(3.72) 

(3.73) 

(3.74) 


we may write G x and G y as 

- a- (- 

• j 8 +j l x ‘. 


(3.75) 


• * 8 v/j 2 + (i-l) 


- J|A 


(3.76) 
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to be used in the actual implementation of the system. Since each of the above 
relations are similar, we are required to store only one of them and alter the signs 
accordingly. It should be noted that, however, this is not the most efficient method 
of storage. Storing only a few of the cross term values and using an interpolation 
scheme will reduce the storage considerably. Of course, an interpolation table of 
(3.75) and (3.76) will lead to a substantial reduction in memory at the expense of 
some computational efficiency. 

Assuming that the positive sign is chosen in equations (3.75) and (3.76), we have 


T+ — n — 
a 2 fri * a 2 

4 4 

T+ — T+ 

a i i>2 ** a i b? 

3 3 

T+ = T+ 

x a 2 f>3 * a 2 &3 

4 4 

T+ = 7 — 

a 1 6 4 a 1 64 

3 3 


Choosing the positive sign for the (3.63) and (3.64), we also find 


s:, t , 

= Ki, 

3 

3 

^<12^2 

1! 

w 1 
w 

4 

4 

c+ 

°a 3 6i 

= S a 3 b , 

3 

3 

°a 4 6 2 



4 4 


(3.77) 


(3.78) 



32 


Thus, the elements of A BI may be written as 


1 +S a 1 b 1 

T+ 

a l f >2 


R x 



T+ 

a 2 fci 

r + s~ h r 0 + 3 



s:, h R, 

T + 

0362 

1 + S * 3 b 3 

T+ 

a 3 6 4 


T.\ 

SI U R y 

T a + 4b 3 


I + S 

0464 


T ai h 

S<nb 3 Rx 

T a 7 l 

u 


T~, 
a 2 &1 

^a 2 b 2 

^0363 

$a 2 b t Ry 


SajbiR-x 

1 

03 62 


7 — 

^0364 


7*4 bl 

$a t b 2 Ry 

^463 




(3.79) 


The elements of the adjoint of A BI required in the implementation of the CG algo- 
rithm are 


(/ + !..)• = 


( L ab D b y = 



(T* b ,r 

«*)■ 

(%,r 

(Tj h Y 


(T^r 




1 +(*:,>,)• 



Rlisur 


i + (s- bt y 


DJ 



&, h Y 

<^,7., )“ 

Rl(s~ b ,Y 

(T-J" 


(T;, b ,Y 

(s* h Y 


K(s:,hY 


R T AS- h Y 

i'b'ajbiY 

(s*>,Y 

07o.fr 3 )“ 


(W 

%(s- it Y 

(W 

(s:. b ,Y 


3.2.4 A CGFFT Implementation 


(3.80) 


The conjugate gradient algorithm presented in section 2.3 is employed for solving 
the FE-BI system (3.55). The required numerical computation of Az and A a : 


z are 
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performed by first decomposing the system into a summation of two matrices; one 
involving operators associated with the boundary integral and another involving the 
elements of the finite element matrix. Then system matrix-vector products may then 
be written as 



( 3 . 81 ) 


( 3 . 82 ) 


( 3 . 83 ) 


( 3 . 84 ) 


( 3 . 85 ) 
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In each case, the operation is performed such that only the resulting vector {s} need 
be stored. 

3.3 Computational Considerations 

The FE-BI method for rectangular enclosures is efficient in terms of memory 
usage and computation time, and each of these aspects is discussed in detail below. 

3.3.1 Storage Efficiency 

The fundamental advantage of this method is the reduction of storage require- 
ments, so that the scattering by electrically large bodies may be evaluated. To show 
that the low storage requirement of 0(N g ) is assured, we refer to Tables 3.3 and 3.4. 
These contain a list of all major memory consuming variables. A summarized list 
is also given in Table 3.5. Specifically, Table 3.5 includes the memory requirements 
pertaining to the finite element matrix (FE), fast Fourier transforms (FT), boundary 
integral cross terms (Cross) and conjugate gradient variables (CG). N c is one less 
than the number of elements connected to a particular node, and a typical value of 
5 is used here. 

To put the quantities of Table 3.5 in terms of N g , the total number of nodes, we 
consider two special geometries. The mesh in Fig. 3.3 corresponds to a penetrable 
body, while that of Fig. 3.4 corresponds to an impenetrable structure tightly enclosed 
by the picture frame. Within each special case two extremes are considered; a mesh 
corresponding to a square object (N x = N y ) and a long strip (N x » N y ). In each 
case, N g is assumed to be large. 

Alluding to Table 3.6 the total storage is 0(N g ) for the square region, but is 
somewhere between 0{N g ) and 0(N J) for the (N x » N y ) case. This is based 
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Memory Consumption 

variable 

type 

count 

comment 

Mesh 


R 

n 3 

X coordinate of global nodes 


R 

N g 

Y coordinate of global nodes 


I 

37V e 

Node-element connectivity 

Pointers 

ABint 

I 

N a 6 

Observation and integration points 

Pnodes 

I 

Pnum 

Nodes on conducting bodies 

dmatl 

I 

N g - N a b 

Element material properties 

Finite Element Matrix (FE) 

Ar 

C 

~ W)(Af. - AT.) 

Non-zero values of FE matrix 

col 

I 

~ (“ipKAf, - N.) 

Column pointer of FE matrix 

row 

I 

N g -N a 

Pointer to rows of FE matrix 

Conjugate Gradient (CG) 

Phiz 

C 

N g 

Unknown vector 

CGI 

c 

N g 

Residual vector 

CG2 

c 

N g 

Search vector 

CG3 

c 

N g 

Temporary vector 

q 

c 

MAX(N x , N y ) 

Temporary vector 

phiinc 

c 

K 

Incident field vector 


Table 3.3: List of major memory-consuming variables 
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Memory Consumption (continued) 

code 




variable 

type 

count 

comment 

Fourier Transforms (FT) 

FTHxl 

c 

to 

Fourier transform along x-direction 

FTHx2 

c 

to 


FTHx3 

c 

2N X 


FTHx4 

c 

2N X 


FTHyl 

c 

2N y 

Fourier transform along y-direction 

FTHy2 

c 

2Ny 


FTHy3 

c 

2 Ny 


FTHy4 

c 

2N y 


FT 

c 

2MAX(A r I , N y ) 

FT of unknown sub-vector 

WR 

R 

2UAX(N x ,N y ) 

Temporary array 

WI 

R 

2MA X(N x ,N y ) 

Temporary array 

Cross-Term Matrices (Cross) 

PQp 

C 

~ MAX(A' x , N y ) 


PQm 

C 

~MA X{N x ,N y ) 



Legend 


R = REAL*4 
C = COMPLEX 
I = INTEGERS 


Table 3.4: List of major memory-consuming variables (continued) 












37 


Major Memory Consumption (N T > N y ) 

Item 

Type 

Count 

FE 

COMPLEX 

(^)[JV. - 2(JV, + A',)] 

FT 

COMPLEX 

127V X + 8N y 

Cross 

COMPLEX 

2 N* 

CG 

COMPLEX 

AN g 


Table 3.5: Summary of major memory consumption 


Major Memory Consumption: Penetrable 

Item 

II 

N x » Ny 

FE 

(^)( 

i *,«) 

FT 

20 ^ 

12 N,/(N, + 2) 

Cross 

2N g 

2( NJN y f 

CG 

4 N g 

4JV S 

total 

~ 9Ng 

EHUUSSH 


Table 3.6: Summary of major memory consumption: filled mesh 
























Figure 3.3: Example of the mesh of a penetrable structure 



Figure 3.4: Example of the mesh of an impenetrable structure 
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on the assumption that all cross terms are individually stored, but by using an 
interpolation table, the 0(N g ) memory requirement can be assured regardless of the 
value of N x with respect to N y . In Table 3.7, more dramatic results for the storage 
of the cross term are listed. In this case, all of the unknowns are on the outer two 
boundaries, so it appears that the storage is 0(N £) for the square case. One must 
note, however, that the factor multiplying the N g term may be quite small. The strip 
case, on the other hand, requires an 0(N g ) storage. This case would be an unlikely 
candidate for the use of this method, since that structure would be handled much 
more efficiently via a direct implementation of the CGFFT method. As noted above, 
the storage of the cross terms may be brought down to 0(N g ) for all cases by using 
an interpolation table, and this will certainly be necessary in a 3-D implementation. 
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Major 

Memory Consumption: Impenetrable 

Item 

ii 

N x » N y 

FE 

(^■)N ,/ 2 

)N,/2 

FFT 

5N g /2 

CO 

£ 

Cross 

N]/32 

n;/s 

CG 

4^ 

4 N g 

total 

~ A/^/32 + 8 N g 

~ N*/8 + l7N g /2 


Table 3.7: Summary of major memory consumption: open mesh 

3.3.2 Computational Efficiency 

Since the primary importance of the FE— BI method is storage reduction, a com- 
parable level of efficiency with alternative methods is a bonus. A method for which 
a fair comparison may be made is the standard CGFFT. This requires a 2-D FFT, 
which is slower than using multiple 1-D FFTs for large bodies. We compared the 
two methods for a specific penetrable scatterer using an Apollo 3500 without code 
optimization. The pertinent CPU times are compared in Table 3.8. The comparison 
provides only a relative measure of the speed difference. 


3.4 Far Field Computation 


The scattered fields may be computed as 


nv) = - j r {G(p,f) 


-m 




dl' 


(3.86) 
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Body Properties 

FE-CGFFT 

CGFFT 

Composition 

Dimensions 

T/I (s) 

■ 



fl 

Total 

dielectric 

2A x 2A 

8.6 


D 



5610 

e r = 4-j.l 



■ 

■ 





Legend 


T/I = time/iteration 
I = number of iterations 


Table 3.8: A comparision of computation efficiency of the FE-CGFFT with the 
CGFFT method 

Using the discretization scheme developed earlier, we have 
<f> a (x,y) = 

~ { Jr [ - K+G(x,x',y,y Cl )<f>(x',y bl j\ dx 
+ L [ K t G ( x ’X C3 ,y,y')(f>{x a2 ,y') - K-G(x,x cj ,y,y')<j>(x h ,y')]dy' 

•'l C2 ■* 

+ L [ K yG(x,x\y,y C3 )<l>(x\y a3 ) - K~ G(x,x' ,y,y Cz )4>(x' ,y b3 )\ dx 1 

c 3 J 

+ J r [ K x G ( x i X C4 ,y,y')4>( x <n , y') - Kt G( x , x <* , y , y')<t>{x hA , y')] <V j 

(3.87) 

where the definitions for and are as specified previously. Letting 

1 f x j + ^ 

/3 r (i,y,y c ) = 7 - / A G(x,x',y,y e )dx' (3.88) 

0 y ( x ,x c ,y) = L f 3 2 G(x,x c ,y,y')dy' 

^ J v 3 - f 


(3.89) 
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7 x (x,y,y c ) = l - 

r^j+y Q 

I Xj .a QjG( x , x ',y,yc)dx' 

(3.90) 

7 v (x, x c , y) = i 

fiO + y d 

4 -» 

(3.91) 


(3.87) becomes 


<f> a (x,y) = 

( Nx 

-|E ([^ r ( x > J6 Vc } ) - 7 r (x, y , j/ ci )] {<&., } j - [0 x {x, y, y Cl ) + 7 x (x, y , y Cl )] {^ 6] } ; ) 

N y 

+ E ([^ V ( x ’ x ^> y) + 7*(x, ar C2 , J/)] {<£a 2 }j - [£ v (x, x C2 , y) - l y (x , x C2 , y)] {<fr, 2 } •) 

i=i 

N x 

+ E([^ I ( x >J/'J/«) + 7 I (a;,y,yc3)]{^a 3 }i -[/3 x (x,y,y C 3)-7 x (x,y,y C3 )]{^ fc3 } J ) 

i=l 

'j 

+ E (^ y ( x ’ X C4 » y) - 7 y (x, X ct , y)] {<^a 4 }, - \P y {x, x Ci , y) + 7 y (x, x C4 , y)] {<f > bi },) 1 

(3.92) 


valid for all observation points (x,y). To specialize (3.92) to far zone computations, 
we must introduce the appropriate asymptotic expansion for the Hankel functions 
implied in (3.88)-(3.91). In doing so, we have 


where 


/T(x,y,y Cl ) = jf 0 (p)M6,y Cl )e jfc ^ c “* 

3 3 

(3.93) 

P v (x,x C3 ,y) = jf o (p)f2{0,x Ci )e jk °y>* iae 

< 4 

(3.94) 

7 x (x,y,y Cl ) = -fo{p)fi(6,y Ci )k o Asm0e ]koX > cose 

3 3 

(3.95) 

7 v (x,x C2 ,y) = —f o {p)fi{0, x C2 )fc 0 A cos 6e’ k ° yi sin * 

< 4 

(3.96) 

- UxKf*" 

f (a x “ nfi 

/i(^,y Cl ) = -e 3 

3 

r /a \ ifcoI ‘ : 2 cosff 

/2(0,x C2 ) = -e 

4 

(3.97) 

(3.98) 

(3.99) 


(3.98) 

(3.99) 
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in which (p,0) imply the usual cylindrical coordinates of the observation point. Sub- 
stituting expressions (3.93)-(3.96) into (3.92), we obtain 

<f>jf(x,y) = -Mp) 

( 

jE(L? + A: 0 A sin 6»] {<f> ai } : - [j - k 0 Asin0] {<^ 6 , };) h{9, y Cl )e jkoT > cos(? 

Ny 

+ E (U ~ k ° A cos Q] {$*2 }] - [j + cos 0\ {<t>b 7 }j) f 2 {0 , Xc 2 )e jkoy > sin 8 
j=i 
N z 

+ E ([j - k 0 As\n9] {Mj -\j + k 0 A sin 6} {^3) M9,y C3 )e jk ° x ’ cose 
i=i 

Ny \ 

+ E ([j + k °A cos 6) {<}, -\j- k a A cos 6} x Ct )e^ sin6 i 


(3.100) 


From (3.100) the echowidth becomes 


a = 


4k a 


N* 


E(L? + k °A sin 6] {4, ai } 3 - [j - k 0 A sin 9} {<f > 6 , },-) /, ( 9 , y ci )e ifc ^ cos<? 

3 - 1 

Ny 

+ E (L? - ^Acosd] - L? + &<> A cos 0] {&,},) / 2 (*, *c 2 )e jfcoS " sin * 
j = l 

+ E(U - *»Asin0] {*.,},• - [7 + A : 0 A sin 0] {^i)/i(*,yc 3 )e i *»*> co “* 

3 = 1 
Ny 

+ E (L? + cos 0} « }, -\j- k 0 A cos 9} {<f> bt } ] )f 2 {0, x Ci )e ]k ° y > sin 

i=i 


(3.101) 


3.5 Code Validation 


The scattering patterns from several rectangular structures are presented. The 
echowidth is computed for each structure and compared to the results of the moment 
method. The bodies are discretized at a sampling rate of 20 samples/free-space 
wavelength. 

Results are presented for the following cases: 
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• perfectly conducting bodies (Figs. 3.5 and 3.6) 

• partially and fully coated perfectly conducting cylinders (Figs. 3.7 - 3.12) 

• composite body (Fig. 3.13) 

In each figure, the discretized geometry is shown along with the corresponding 
results. As seen in all cases, the generated patterns using the FE-CGFFT formulation 
are agree with the moment method data. 
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Backscatter Echo WidthA 
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.25 x 1. X Coated Perfect Conductor (E-pol) 



Angle (deg) 


Figure 3.7: E z backscatter from a .25 x 1A perfect conductor with a A/20 thick 
material coating containing the properties t T = 5. — j. 5, /i r = 1 
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.25 x 1 . X Coated Perfect Conductor (H-pol) 



Angle (deg) 


Figure 3.8: H z backscatter from a .25 x 1A perfect conductor with a A/20 thick 
material coating containing the properties e r = 5. - j. 5, /x r = 1 
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m 


m 

*SM 



25 x 1 . X Coated Perfect Conductor (E-pol) 
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.25 x 1. A, Coated Perfect Conductor (H-pol) 


FE-GCFFT 



0.0 15.0 30.0 45.0 60.0 75.0 90.0 


Angle (deg) 


Figure 3.10: H z backscatter from a .25 x 1A perfect conductor with a A/20 thick 
material coating containing the properties t T = 5. - j.5,/i r = 1.5 - j . 5 



Backscatter Echo WidthA (dB) 
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.25 x 1. A, Double Topped Conductor (E-pol) 

30.0 

20.0 
10.0 
0.0 

- 10.0 
- 20.0 
- 30.0 

- 90.0 - 60.0 - 30.0 0.0 30.0 60.0 90.0 

Angle (deg) 



Figure 3.11: E z backscatter from a .25 x 1A perfect conductor with two A/20 thick top 
material coatings. The lower layer has the properties t r = 2. — j. 5, fi T = 
1.5 -j.2, and the upper layer has properties e T = 4. - j.5,fi r = 1.5 -j. 2 
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.25 x 1. X Double Topped Conductor (H-pol) 



Angle (deg) 


Figure 3.12: H z backscatter from a .25x 1A perfect conductor with two A/20 thick top 
materia] coatings. The lower layer has the properties e r = 2. — j. 5, = 

1.5 — j. 2, and the upper layer has properties t r = 4. — j. 5, /x r = 1.5 — j . 2 




.5 x 1 . X Composite Body (H-pol) 



0.0 30.0 60.0 90.0 120.0 150.0 180.0 


Angle (deg) 


Figure 3.13: H z scattering from a composite body. Both the perfect conductor and 
the dielectric body are A/2 in each dimension. The material properties 
are t T = 5. — j. 5, p T = 1.5 — j . 5 
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3.6 Summary 

A procedure was developed for computing the scattering by 2-D structures. This 
procedure combined the boundary integral and finite element methods, and the re- 
sulting system was solved via CGFFT. The main advantage of the new approach was 
a reduction in memory demand to 0(N) compared to 0(N 2 ) required with tradi- 
tional solution techniques. A detailed map of the storage requirements was presented, 
and the principle memory consuming arrays were discussed. Also, the computational 
efficiency of the technique was examined and found favorable. To validate the pro- 
posed solution approach, several backscatter patterns were presented and compared 
with results based on traditional solution methods. 



CHAPTER IV 


A Finite Element — Boundary Integral Method 
for Two-dimensional Scattering Using Circular 
and Ogival Termination Boundaries 

4.1 Introduction 

It is possible to choose other boundaries that result in convolutional integrals, 
and in this chapter we consider circular and ogival enclosures. Clearly, a circular 
enclosure would be attractive for circular scatterers whereas an ogival boundary will 
be more attractive for those structures conforming to this boundary. In the case of 
the circular boundary the entire integral is convolutional ensuring the O(N) memory 
demand of the system provided an iterative solver is used. When an ogival enclosure 
is used the integral becomes convolutional only if the observation and source points 
are on the same arc, but an efficient storage scheme is again required for the remaining 
“cross-terms” 1 . 

In the following sections, the pertinent FE-BI formulations are developed for the 
circular and ogival boundaries. Results for several circular and ogival structures are 
presented and shown to be in excellent agreement with that obtained by traditional 
methods. 

1 “cross terms” refer to integrals for which the source and observation points are not on the same 
arc 
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y 



Figure 4.1: Geometry of the scatterer 

4.2 Analysis 

Consider the plane wave 

7 nC (p) = z<f>T c (p) = ze jk ° pcoa l e - e ° > (4.1) 

illuminating a composite cylinder as shown in Fig. 4.1 and we are interested in 
computing the scattered field. For the application of the Finite Element - Boundary 
Element Method the target is enclosed in a fictitious circular or ogival boundary as 
shown in Figs. 4.2 and 4.3. Within the boundary r a , the finite element method is 
used to solve the Helmholtz equation 


V • + klv(p)<f>(p) = 0 (4.2) 


where 


M/>) 


4>{p) = E 2 (p), u(p) = 


v(p) = t r (p) 


(4.3) 
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for E-polarization and 

<f>{p) = H z (p), u(p)=-^—, v(p) = p T (p) (4.4) 

trip) 

for H-polarization. The free-space wave number is k 0 = L 0 y/fi o t 0 and fi r and c T are 
the relative permeablility and permittivity, respectively. On the boundary T a the 
Helmholtz integral equation 

m = *~(» - 1 -#?.) }<«. (<-s) 

provides the required boundary constraint, implicitly satisfying the radiation condi- 
tion. In (4.5) 

G(?,p.) =—«?’(*.!? -?.l) (4.6) 

is the 2-D free space Green’s function where H^(-) denotes the zeroth order Hankel 
function of the second kind. Also, denotes differentiation with respect to the 
outward normal, whereas ~p and ~p a are the usual source and observation points, 
respectively, and 

I P~Pa\ = \/(^-^a) 2 + (j/-J/a) 2 (4-7) 

4.2.1 Circular Enclosure 

Discretization of the Scatterer and Field Quantities 

The region enclosed by r a , denoted as R a , is discretized into N e finite elements 
as illustrated in Fig. 4.2. In the figure, p a is the radius of the circle and a a is the 
integration angle along this boundary (Further definitions for the finite element mesh 
are indicated in Table 4.1, while the definitions of the field vectors are indicated in 
Table 4.2.). We note that nodes along r a are equispaced with angular displacement 
A. 
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Table 4.1: 



Figure 4.2: Partially discretized body in a circular enclosure 


Symbol 

Description 

N n 

number of nodes in the finite element mesh 

N g 

number of unknowns 

N e 

number of elements in the finite element mesh 

N a 

number of nodes or elements on T a 

N d 

number of nodes or elements on Vj 


Definition of various quantities associated with the finite element mesh 
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Symbol 

Description 

<f>a 

field at the nodes on T a 

‘Ipa 

normal derivative of the field at the nodes on T a 

h 

field at the nodes region I enclosed by T a and Yd 

<t>d 

field at the nodes on p 


Table 4.2: Definition of various field vectors associated with the finite element mesh 
and its boundary 

Derivation of the Finite Element Matrix 


The weighted residual expression over each element may be written as [5] 

jj R'W'dS' = 0 .' = 1,2,3 (4.8) 


where 


u ( x >y)fc;<f>'{x,y) 


d_ 

dy 


u{x,y)-^4>% x ,y) 


~ k 2 0 v(x,y)4> e (x,y) (4.9) 


In (4.9), W t e is the ith weighting function associated with the eth element. Substi- 
tuting (4.9) into (4.8) and invoking the divergence theorem yields 

+ tfvpW') dYt e 



—u 


dtf_ dWJ_ dtfdWf 
[ dx dx dy dy 


+ 


[ w t 'ip e dr e = o 

Jr* 


where P denotes the contour enclosing the eth element. Additionally, 

M' 


%l) e = u e 


dn 


(4.10) 


(4.11) 


is zero outside element e. Summing over N e elements we obtain 


N c 

E 

"'S ! 



—u 


dp dWf d<f> e dW t e 


[ dx dx dy dy j 

N a . N d 


+ k*v4,’W'\ 


rffi‘ 


'’o r I'd r 

+ £ L w;rdr a + £ / w;rdr’ d = o (4.12) 

*=i Jl ° »= 1 Jr d 
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where the summations over s refer to the elements with sides adjacent to the fictitious 
(r«) and conducting (Tj) boundaries. 

The integral over the conducting boundary in (4.12) vanishes all cases. This is 
obvious when no conductor is present. When <fr = H z , the normal derivative of the 
field is zero on the conductor and the field unknowns on the boundary are allowed to 
“float” (i.e., they are not fixed to a specific value but are determined by solving the 
system). Finally, when <f> = E z , imposing the Dirichlet condition after assembling the 
finite element system results in the elimination of those equations associated with 
the integral over Tj. 

Proceeding with the discretization, the field and its derivative within each element 
may be expanded into a linear combination of shape functions 

? = (4.i3) 

j = 1 

= EWq (4.14) 

Jt=l 

Substituting (4.13) and (4.14) into (4.10) and choosing Wf = Nf (Galerkin’s method), 
we obtain 


where 




j=l k=lj=l 


(4.15) 


and 



dN, e dN '? 
i L 4. 

dx dx 


dN'dNf 
dy dy 


klvNfN 


>} 


dn e 


(4.16) 


J'N'N'dl e (4.17) 

For linear triangular elements, Nf are given by 

(Vf = + 6*x + <y) 


(4.18) 
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with 


W m 

(4.19) 

a l = x)vl - x'kVj 

(4.20) 

b i = y e j - vt 

(4.21) 

c i = x k “ x ) 

(4.22) 


and (x',J/f) being the coordinates of the ith node of the eth element. From (4.18) 


dNf __ b\ 

dx ~ 2 Cl' 
dN[ _ _c[_ 

dy 2fl e 


(4.23) 

(4.24) 


Using (4.23), (4.24) and the identity 


jj ( N;y(N' 2 yjxj y 

S e 


= 2 n e 


( p + <7 + 2)1 


in (4.16) reduces to 


(4.25) 


where 


“h = + ^) - kiv^a + s,,) 


(4.26) 


6 


v 


1 if i = j 


(4.27) 


^ 0 otherwise 

We note that in deriving (4.26) we have assumed that u and v (the reciprocal of the 
material constitutive parameters) are constant within each element and are given by 
u e and v e , respectively. 

To find an algebraic expression for b 3 ik , we may reparametize the integral in (4.17) 


as 


b 


■,=r 

Jot 1 


P’P a k r a da 


(4.28) 
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where P’ and P 2 * are given by 


P i( Q ) = 1 - “I 


a — a; 


oft - a 


Ptt«) = -7 


a — a; 


Qo Qi 


Integrating, we have 


b’k = + 1) 


(4.29) 

(4.30) 

(4.31) 


(4.32) 


Substituting the previous equations into (4.12) a sparse matrix is obtained for 
the nodal fields that has the form 

A aa A a j 0 —B a 

Ai a An An 0 

0 Adi Add 0 

0 0 0 0 

In this, the values of the elements in the submatrix A pq are the contributions asso- 
ciated with the nodes in group (region or boundary) p which are connected directly 
to the nodes in group q. Also, 





0 


<fri 


0 


<frd 


0 


(fra 


0 


[■®aa]ifc — ^2 b* k — — — + 46,* + 6i+i,fc) (4.33) 

«=1 D 

The last row in (4.32) has been intentionally left blank to imply a need for another 
set of equations relating the fields and its derivatives on r a . This additional set of 
equations is produced by discretizing the boundary integral equation. 


Evaluation of the Boundary Integral 

The boundary integral in (4.5) may be rewritten in cylindrical coordinates via 
the transformations 

\ p - pa \ = |*e(p cos Of p a cos cn a ) y(p sin o p a sin ot a ) | 

= + 2 />/>„ cos(a - a,) 


(4.34) 
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where ( p,a ) and (p a ,ot a ) are the usual source and observation points in cylindrical 
coordinates. For \p\ = |p a |, 


1 9 ~ Pa I = 2p|sin(^)| 


(4.35) 


and the Green’s function and its normal derivative may be written as 

G(p,p a ) = - J -H { 0 2) (2k 0 p a sm {?=?*)) 

= ^#j 2) (2fc 0 p a sin(^))sin(^). 

We may now write (4.5) as 

\4>{p,ct) = 4>' nc (p, a) - / 0 (p,a) + f\{p, a) 
where as a result of (4.36) and (4.37) 

1 f 2 tt 

fo{p,a) = -~Pa rpipa,a a )H^ ] (2k 0 p a sin C 2 ^)) da a 

j r 2tt 

fi(p,ot) = -^pa j- <j)(pa,a a )Hi (2k 0 p a sm{2^*))sin {*=?*) da a 


with 


Tpipa,C* o) — ^ &iPa,&a) 

Op a 


(4.36) 

(4.37) 


(4.38) 


(4.39) 

(4.40) 


(4.41) 


The factor of | in (4.38) accounts for the singularity associated with and the 

■f (4.40) denotes principal value. 

We may now discretize (4.39) by expanding the field using pulse basis functions 


as 


N a 


where 


xl>(p a ,a a ) ~ ^P A («a - 

J=1 


1 |Q Q — aj|<Y 

p a K -<*,•) = ; 2 

0 otherwise 


(4.42) 


(4.43) 
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and A is the angular width of the integration cell as indicated in Fig. 4.2. Thus, the 
discrete version of (4.39) may be written as 

fo{p,a) = Yli’j [ 2 (2k 0 p a sin ( 2 ^)) da a (4.44) 

4 j = 1 Ja l~ 2 

Performing point collocation and letting u' — a — a„, we have 

U(P, a,) = (4.45) 

^ j = l 

which may be written in compact form as 


jo ' v ° 

fo(p, a,-) = — j- ]T - ctj) 

4 j=i 


(4.46) 


where 


M a » - a j) = f * H^ 2) (2 k 0 p a sin (f )) du 


(4.47) 


It is clear that (4.46) is in the form of a discrete convolution and can thus be written 


fo{p,a) = DFT- 1 {DFTty) • DFT{h 0 )} 


(4.48) 


where the elements of ho are given by 


h 0 {pA) 


_ f A{l-;|[l„(^)-l]} p = o 


Cu! ™ ( 7 )) <*“' P = 1 A'. - 1 


(4.49) 


where in developing the p — 0 term, the small argument approximation of the Hankel 
function was used and is given by [12] 


lim #«>(*,,)= I 1 ^ " ° 

P &T- + , 2"(n-l)! j 

t 2"n! T J ■wlkp )" n > L 


(4.50) 
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in which 7 ~ 1.781. Through a similar analysis, the field may be approximated by 
the expansion 


N a 

<f>(p a ,Q a ) ~]T]P A (a a - Q j)4>j (4.51) 

i=i 

and by substituting this into (4.40), we obtain 


i. n N a 

flip, a,) = J—j 1 <t>jhi(a!i - Qj) 


j=i 


(4.52) 


where 


> , , r(°>- Q j )+ f m 

h\(a,-Q J )= l J (2/r 0 /? a sin (^)) sin (J^) 

J (o.-oj)-f 

Clearly, (4.52) may again be written in operator form as 


(4.53) 


flip, a) = DFT- 1 { DFT(<f> ) . DFT{h x )} 


(4.54) 


where 


hi(pA) = 


ko Pa (f - sin f) +j 


ir^OPa 


P = 0 


( Cljl (2 k oPa sin (£)) sin {%) du> p = 1, ..., N a - 1 
where again (4.50) was employed. 

Point matching (4.38) at each node results in the system 


(4.55) 


\4>i = <t> i r c -f 0 (p,a i ) + f 1 (p,a i ) 


(4.56) 


which may be written in operator form as 


where 


^aa^a T aa 0 a — (f> a 


[Laa\ij — ^ Mat “ aj) 


(4.57) 


4 


(4.58) 

(4.59) 
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Figure 4.3: Partially discretized body with an ogival enclosure 
A final system is obtained by combining (4.57) with (4.32) to yield 

A*a A*l 0 B aa (f) a 0 

Ai a An Au 0 <j>i 0 

= (4.60) 

0 Adi Add 0 <j>d 0 

M aa 0 0 -L aa J [ J C c 

which can be solved via the conjugate gradient algorithm to obtain the nodal fields. 

4.2.2 Ogival Enclosure 

Discretization of the Scatterer and Field Quantities 

The region within r„, denoted i? Q , is discretized into N e finite elements and a 
partial discretization is shown in Fig. 4.3 for the circular case. With respect to 
Fig. 4.3, the definitions of the various quantities are found in Table 4.3. Further 
definitions for the finite element mesh are indicated in Table 4.4, and the field vector 
definitions are indicated in Table 4.5. 
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Symbol 

Description 

A P 

angular displacement between nodes on r ap 

Pap 

radii of T 0p 

Q °P 

angular integration variable along T ap 

t 

distance between centers of curvature of T„ 

a P 

Vc p 

y-coordinate of the center of curvature of T n 

Up 


Table 4.3: Geometric quantities in reference to the figure above 


Symbol 

Description 

N n 

number of nodes in the finite element mesh 

N g 

number of unknowns 

K 

number of elements in the finite element mesh 

N a 

number of nodes (=N al + N a2 ) on T a 

r a 

r -1- r 

Mi T 


Table 4.4: Definitions for the finite element mesh with an ogival enclosure 


Symbol 

Description 

<t>a p 

fields corresponding to the nodes on T ap , p= 1,2 

^a p 

fields corresponding to the midpoints of the nodes on T ap 


fields at the nodal midpoints on T ap 

<t>l 

fields corresponding to region I enclosed by T a and T* 

4>d 

fields corresponding to the nodes on the T d 


Table 4.5: Definition of the field vectors for the ogival boundary 
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Derivation of the Finite Element Matrix 

The derivation of the finite element matrix follows that described in section 4.2.1 
with the exception of the matrix B aa . Consider the ogival boundary as indicated in 
Fig. 4.3. The boundary contour r a is comprised of two arcs labeled r a , and r aj , 
which form the vertices of the ogive where they meet. At the vertices the unknown 
normal field is discontinuous and will therefore be evaluated at the midpoint. Also, 
in evaluating the contour integral, the field derivative are expanded in terms of pulse 
basis functions, rather than linear functions. This results in a different B aa matrix 
and involves the replacement of P* in (4.28) by the pulse basis function expansion 

[ 1 ifO < |a — Qj| < f 

P^(Q-aj)=j (4.61) 

I 0 otherwise 

By integrating in cylindrical coordinates we then obtain 

^ = ^ + *..j+i), 3 = 1, *' = 1,2 (4.62) 

where l e is the length of the eth boundary element along r„ and is equal to p ap A p 
for r 0p , p = 1,2. Performing a summation over all boundary elements then yields 

[Baa\ij = = + ^«J+l) (4.63) 

e=l Z 

where V is the length of the jth element since the jth “node” (associated with the 
unknown iftj) is at the center of the jth boundary element. 

The remainder of finite element analysis for this case proceeds exactly as in section 
4.2.1. 

Evaluation of the Boundary Integral 

The evaluation of the boundary integral along an ogival contour is similar to 
that described for the circular boundary. For integration and observation points 
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on the same arc, the integrals become convolutions. On the other hand, when the 
integration and observation points reside on different contours, the integrals have no 
special form and must be discretized and stored in memory as efficiently as possible. 

The distance between the source and observation points in terms of cylindrical 
coordinates for points on the same arc is given by 

\p - Pa p \ = \fp* + H v - 2 pp a „ cos(a - Q ap ) p=l,2 (4.64) 

When the source and observation points are along different arcs, (4.64) becomes 
I P <1 - Pa p I = \J (p cos <*9 - pa COS a a „) 2 + (p sin a, - ~p a sin Q ap + y c? - y Cp ) 2 

p,q= 1,2 (4.65) 

in which the subscript a p refers to the integration coordinates along contour p and 
the subscript q refers to the observation coordinates. Also, y Cp is the y-coordinate of 
the center of curvature for contour p for p = 1,2. For further reference we note that 
(4.65) may be also rewritten as 


I Pi ~ Pa 2 1 = 

1 

J(p* +Pl 2 ~ 2 PlP °2 cos(oi — ct a2 ) + f 2 2t(pi sino. - p a2 sina 02 ) (4.66) 
V ^ 1 2 1 ^ 1 22] j 

where t = y Ci — y Ci and the upper sign corresponds to the upper set of subscripts. 

To discretize (4.5), the fields are expanded as 



N ai 

N a 





■) P &( Q a - <*j 

Wj+h + £ 

PA(<X a ~ 

" Q Mj+1 

(4.67) 


J = 1 

* 

Vj. 

II 

+ 





N ai 

N a 




HPa, 

a Q ) ^ ]T p a ( a a ~ 

+ £ 

PA(ot a - 

- Ctjtyj 

(4.68) 


i=i 

j-N ai +1 




where as before 








r 

1 if |q q — a, 

\<T 




£ 

P 

p 

1 

P 

«-». 

ii 

1 U J | 

1 — 2 


(4.69) 



0 otherwise 
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and 


^i+i — + ^J+l) 


(4.70) 


Substituting (4.67), (4.68) and (4.69) into (4.5) then yields 




“ / G 0 (/»i,/» Ol ,Oi -« a] )p ai da ai 

j=i 7o j 


+ 


01 rofj + Ai ^ 


Wo, 

j=i 7 «> 


dp ai 

£>j+A 2 


G 0 (^>1 , p ai , Qx - Q a , )p a , da ai 


^ raj+*2 

2 ^ i I G 0 [p\ , Paj 7 ®1 1 ^02 )^aj doia^ 

j=Na, +1 7 °J 

^2, r a ]+ A 2 3 

d” / ^ ^j+r f a ^o(/^l i /^a2 ) ^1 j )/^aj 

i=N«,+l 7 ° 7 >a 2 

when the observation point is on r a , and 


(4.71) 


2^(/ , 2,<>2) — 4>' nC {p2,C*2) 


I G 0 [p2y Pa \ 5 ^ 2 ) ®oi )^oi do; a , 
>=1 7o 0 


+ 


“• rQj+Ai Q 

^j+h T o (*o(P2) Pa\ i &2i Q; a] )Pai d<J 0 ] 

j=l 2 7 "J 


dp ai 

raj+& 2 


raj+Ai 

~ 2_rf / G 0 (p2,Pa 2 ,Ol2 ~ a a2 )pa 2 dO‘a 2 

]=Na x +1 70r J 

ra,+A2 Q 

+ Z-, ^j+h T G o{P2,Pa 2 ,Ot2 ~ Oa 2 )pa 2 da a2 (4.72) 

i=Noj+l 7 a J 

for observation on r 02 . Performing point collocation at the nodal midpoints, (4.71) 
and (4.72) further reduce to 


!#/*.«* +*) = ^(Pi.Oi+j) 

N. 


y=l /*Oi-a; + f 

Go{pi,p ai ,u)p ai du 

j = 1 ■'o.-orj + f 
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; ra t — Qtj + y 

+ z2^j+hr — G o(pup ai ,u)p ai du 

j=l 2 j o.-Oj + f Op ai 
^ fOj+^2 

x > / ^ol/^l -> Pa? •> ^i+ 1 1 ^02 )/^Q2 doi a2 

J=N„ 1+ 1 •'“J 

rcifj+Aj ^ 

+ zl 4>J+). f -^—G 0 {p u p a2 ,a i+ x a a2 )p a2 da a2 (4.73) 

3=N ai +1 J °) °Pa 2 2 


for observation on r oi and 


\4>{p2, a ,+ i) = (f> ,nc (p2,a l+ L) 


r^+Ai 

/ t ^ j I G o(P2i Pa\ i > ®aj ) Pa] dck ai 

j=l •'“J 2 

No, 

2 /■a, 

+ em/ 

i=l J 


a_,+Ai ^ 


Pai i ^,-+1, ^O] )/^ai do aj 


rai 

ra,-a ; + £- 

- IJ V’i / . Go(p2,p a2 ,u)p a2 du 

}=N ai + 1 -'Ori-Oj + y 

^ ra.-otj + f Q 

+ ^j+J T .a — G 0 (p 2 ,p a2 ,u)p a2 du 

j=N ai + 1 7 v.-Cj + fOp^ 


(4.74) 


where the in the subscript refers to the fictitious “node” midway between the 
actual nodes. A system of equations can now be obtained by testing (4.73) and 
(4.74) at a sequence of points on the contours. This yields 


4 G\ 4> ai — — |Th0 Oi -f P\ 2 l}l a2 — \M\\C\<}> ai — ^Ql2C 2 <f> a2 } 

\C 2 <j>a 2 = <t>'a 2 + L ~ { ^21 + ^22^02 “ — 2 ^22^2^2} 

which can be alternatively written as 


-M11C1 

Q21C1 


Q12C2 

M22G2 






- 


" 




<£aj 



£11 

P12 


4 > a 2 


A ,nc , 

^“i + j 

P21 

L22 




^“ 2 + J 




1 

O 

W 

1 




(4.75) 


(4.76) 
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In this, the nonzero values of the upper bidiagonal matrix C are 1, and 4>' nc , the 

°P + 2 

value of the incident field evaluated at the nodal midpoints. The matrix D accounts 
for the double use of the nodes at the endpoints and the remaining elements are 
given by 


•Mpp — |(P ~ Mpp) 

Qpq = — ?Qpq 


for p — 1, 2 in which 


(4.77) 

(4.78) 

(4.79) 


= 

/■a,-aj + ~ E 


/ i- Go{ppj Pa.pi u )Pa p du 

(4.80) 

(Mi, « 



T Ap a G 0 {pp, Pa p iU)p a du 

J a.-aj + ^f Opa p 

(4.81) 

[Qpdij “ 

rOj+Aq 


/ G 0 (/) p ,p a? ,a t+ i,Q a Jp a 9 dQ a4 

^ Of j • 

(4.82) 


Qp Go{ppi pa q » + i> G ^a q )pa <t ^ Q a q 

More explicitly, upon evaluation of the integrals 

l m Pr>Uj ~ S , , 

CZ% H\ 2) {2k oPap sin |) sin \du i ± j 

‘ 1 u j • 2 

... )-;^{l-;l[ln(^)-l]} i = j 

l L PPUj ~ \ A. 

C-Z+4 sin *)*. i # j 

[ 0 1 = jffWgT) 

l 21J «■» J A L (L./r=TT\ 


[ Pa 2 ro,+A 2 

p \ ?] 0 = 4 1 


(^oV^rT) 

Pa 2 - P\ cos(ai - a a2 ) ± t sin a a2 da a 
1 2 2 1 ij 1 


(4.83) 


(4.84) 


(4.85) 
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where the upper sign corresponds to the upper set of subscripts on P or Q, while the 
lower sign corresponds to the lower set of subscripts. Introducing the definitions 

M\\C\ Q12C2 

Ki = 

Q21C1 M22C2 

Lu P\ 2 

P21 L22 

the system (4.76) may be combined with that derived via the finite element method 
to obtain 

A a a A a I 0 Baa 4*a 0 

A la All Aid 0 0 

= (4.90) 

0 Adi Add 0 <f>d 0 

A'i 0 0 I< 2 rpa 4>Z i 

J L J L J 

We note that (4.90) can be solved via the CG algorithm to take advantage of the con- 
volution operators M and L in reducing the memory requirements. This algorithm 
is given next. 

4.3 A CGFFT Implementation 

In a manner similar to the previous chapter, the matrix vector multiplications Az 
and A a z is represented as a summation of matrices, one corresponding to the finite 
element portion of the system and the other corresponding to the boundary integral 
portion. Thus, a typical product may be represented as 


(4.89) 




{s} = { s } b / + { s}fe 


( 4 . 91 ) 
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where 


Wsi = 


0 0 0 0 
0 0 0 0 
0 0 0 0 


Z\ 

z 2 

z 3 


and 


{•s}fE = 


For the adjoint operations, we have 


[ K x 0 

0 I< 

•n 


A aa A a I 

0 

B aa 



Al a An 

Am 

0 


*2 

0 

0 

Add 

0 



0 

0 

0 

0 


*4 


{•sJb/ = 


*5 

O 

O 

o 


z\ 

0 0 0 0 


z 2 

0 0 0 0 


z 3 

S? 

o 

o 

o 

J 


Z 4 


and 


{•s}fe = 


Aaa A a ar 0 0 


Z \ 

A a ia A)j A) d 0 


z 2 

O 

O 


Z 3 

[ B aa 0 0 0 


Z 4 


(4.92) 


(4.93) 


(4.94) 


(4.95) 


Again, the operation is performed such that only the resulting vector {s} needs to 
be stored. 


4.4 Scattered Field Computation 

In this section the expressions for the scattered field and radar cross-section are 
developed for both the circular and ogival boundaries. The scattered field is com- 
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puted from the identity 




d 


dn. 


4>{p a ) 


<f>{Pa) 


dn. 


-G(p,p a )\ 


dT„ 


(4.96) 


from which the echowidth is calculated. 


4.4.1 Circular Boundary 


The scattered field expression (4.96) may be written as 


4>\p,&) = -fo(p,a) + fi{p, a) 


(4.97) 


where 


fo{p,a) = ~^Pa J o ip(pa,a a )Hf i) (k Q \fp 2 + p 2 a - ‘IpPa cos (a - Q n )j da a (4.98) 


and 


j [ 2r (* 

/i(p,a) = -p a k 0 / <t>{p a ,a a ) — f 

4 J ° \J, 


//J 2) {h\Jp 2 + pl~ 2 pp a cos(o - a„)) 


P 2 + Pi - 2ppa COS (a - Q a ) 

[p a - pcos(ot - Qr a )] da a 


(4.99) 


To evaluate the integrals in (4.98) and (4.99) we invoke the field expansions (4.42) 
and (4.51). We have 


fo{p,a) - -i- Pa I ^ H (2) (koifp 2 + ~p\ - 2pp a cos (a - q„)) da a (4. 

‘i J=1 \ J 


100 ) 


and 


r< ^ n + * H i (* 

U ( P , = T kop a 2^<f > J / 7 

J = 1 Jq J~T yf, 


j+f # 1 (2) [ko\jp 2 + pl~ 2(>Pa cos (a - a 0 )) 


P 2 + Pa ~ 2pp a cos(a - o a ) 

[p a - pcos(a - a a )] da a (4.101) 


where the remaining integrals over the subsections must be evaluated numerically 
for arbitrary observation. However, for far-field computations (p — ► oo), the Hankel 
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functions may be approximated as 

~ M rc ~ ,ip 

and since 


(4.102) 


for amplitude 


\/p 2 + Pa ~ 2 ~pp a C0S(Q - Q a ) ~ { ' (4.103) 

p — p a cos(a — a a ) for phase terms 


(4.100) and (4.101) become 


fo{p, ot) = 


(4.104) 


and 


/l(p,a) = ^^\hM- e ~ jk ° P 22 COS ( O' - Qj ) e ^0P a cos(a- aj ) 


irk 0 p 


(4.105) 


j=i 


Substituting (4.104) and (4.105) into (4.97) we obtain 


,a)= 


4 V nkop 


Na 


iHV’ j e jkopaCO * {a - Q > ) 

L j = i 

Na 

+ko 22 & cos (a - aj)e jk ° Pa cos (“~^) 
i=i 


(4.106) 


and from (2.31) the echowidth is given by yields the echowidth 

(PaA) 2 


a = 


87T 


Na Na 

j J2 ip je jk OP. «*(«-«>) + g 4>j cos (a - a J -)e J "* 0Po cos(o “^ 

j=i i=i 


(4.107) 


4.4.2 Ogival Boundary 

Following the same discretization scheme used in Section 4.2.2, (4.96) may be 
written as 

Na 


f N a 

<t>’{p,<x) = -{ 22 V’.Mp, «,<*;)+ 22 1 Pfn(p,a,atj ) 

i= i i=N a ,+i 

^“1 Afa . 

- 22 £/« (p> «j) - 22 ^-MP, O', aj ) \ 

i = 1 i=N aj +1 J 


(4.108) 
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where 


flpiP) Qj ) — 

[°j+&p 

1 G a (p, p ap , Q , 0 Q p )Popdo Qp 

Jol, 

(4.109) 


J 

rotj + A p Q 

1 o (/^ /^ap i 

■ / a r j C'Pap 

(4.110) 



(4.111) 


in which 


Pap 5 ^ i G-tip) 

~ 4 H ° 2) { k °\fp 2 + flip “ 2 M*p cos ( Q “ S) + yf P ~ 2j/ Cp (/>sin a - p ap sin a Qp )) 


(4.112) 


and 


5 


<9/?a 


~G 0 (Pi Pap i ^a p ) 


_ jfcp ^1 (2) (* 0 j/jj[ + Pap - 2 H°„ cos (<* - Q ap) + ^Cp - 2y Cp (psm o- - p ap sin Q 0p )) 
4 0> 2 + /»5p “ 2 PP* ~p cos ( a - Q ap) + y? p - 2j/ Cp (/)sin o - p ap sin a Qp ) 

[pa p ~ pcos(a - a Qp ) + y Cp sin a„, 


(4.113) 

and Vc p are the corresponding ^-coordinates of the arc T ap . Using the large argument 
approximation for the Hankel function and the approximation 

\fp + P\p - 2 ~pp7 p cos (a - Q ap ) + yf p + 2t/ Cp (psina - p ap sino ap ) 

p for amplitude terms 

P ~ Pa p cos(q — a Qp ) + y Cp sin a for phase terms 
for p — ♦ oo, the Hankel function simplifies to 


(4.114) 






nk o p 


jkopg-jko -pa p cos(a-a ap )-y Cp sin a] 


(4.115) 
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Similarly, 


#l 2) (/iv)) 


M ~ -jJ-^e- jk ° p e- jk °l- pa ’> cos l a - a ^- v ‘p s H cos(q - a ap ) 


ypV^) J \l7rk 0 p 


(4.116) 


Substutiting these into (4.109) and (4.110) and performing midpoint integration 


yields 


fipipl <*5 ) — 


h P (p,a, Qj ) - k 0 4 “ P ^/ 7r ^< 


Thus, from (4.108) 


j ApPap y 2 j — ^ 2 ^) — ycp sin aj ^ 

k ^pP a P J ^ c -jk o 0 c -j k ° [~^ Q p c °s(Q-Qj- )-Vc p sin aj 


/ . \ 

cos K + - Qap ) 


(4.118) 


r *•> r 

{j ^Z^Le-M 


—Pa j cos(o— Qj — j 1 -)— Vcj sin a 


+ jAlpa 2 ]T ^ + ie ‘^°[' Pa2 sin or] 


j=N 0l +1 


01 A 

+ k 0 A lPat £ <f> J+ ie~ 3ko ^ ^( Q - a} -^)-y Cl sina] cos(a _ a . _ 

3 - 1 * 2 

N a . 

+ Ar Q A 2 p a2 2 ^ + ie" j ' Co l" P ° 2 c “(“-^-¥)-Vc 2 sina] cos(a _ aj _ ^2)1 
i=N 0 ,+i 2 2 J 

(4.119) 


the echowidth becomes 


l *-« 

<7 = — jA,p a . V) . J e _ -’ At> r pa i cosla-aj.-J-l-Vc, sino] 

07T Ji "2 

J = 1 

Na 

-j- jA 2 Pa 2 Y1 ^ + ie" 5><, t‘ Paj Sino] 

J=Afa,+l 2 

N., 

+ k 0 A } p ai Yi <f> j+ Le~ jko l~ p ^ ”*°- a 3-¥)-*i S H cos^ + — - 
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N a 


+ k 0 A 2 Pa 2 £ 008 ( 0 , + — — o aj ) 

j=w 0 ,+ i 2 2 


(4.120) 


4.5 Results 


The scattering patterns of several circular and ogival cylinders for both E- and 
H- polarization are shown in the figures to follow. Figs. 4. 4-4. 6 contain circular ge- 
ometries both coated and uncoated, while Figs. 4. 7-4. 9 contain coated and uncoated 
ogival structures. The echowidth is computed for each structure and compared to 
the results of the series solution for the circular bodies and moment method [13, 14] 
for the ogival structures. As seen in all cases, the generated patterns via the FE-BI 
formulation are in excellent agreement with the corresponding data based on the Mie 
Series and Moment Method Solutions. 
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Perfectly Conducting Cylinder R=.5 X (E-pol) 



Perfectly Conducting Cylinder R= 5 X (H-pol) 



Figure 4.4: E z and H z bistatic echowidth of a perfectly conducting circular cylinder 
of radius 0.5A 
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Coated Conducting Cylinder R=.5 X (E-pol) 



Coated Conducting Cylinder R=.5 X (H-pol) 



Angle [deg] 


Figure 4.5: E z and H z bistatic echowidth of a perfectly conducting circular cylinder 
with a conductor radius of .5A and a coating thickness of .05A containing 
material properties t T = 5 — j 5, fi T = 1.5 — j’0.5 



82 


Coated Conducting Cylinder R=3 \ (E-pol) 



Coated Conducting Cylinder R=3 A (H-pol) 



Angle [deg] 


Figure 4.6: E z and H z bistatic echowidth of a coated circular cylinder with a conduc- 
tor radius of 3A and coating thickness of 0.05A with material properties 
£r 6 — 1.6 j0.5 
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.5 x IX Coated Conducting Ogive (E-pol) 



.5 x IX Coated Conducting Ogive (H-pol) 



Angle [deg] 


Figure 4.8: E z and H z backscatter echowidth of a .5 x 1A perfectly conducting ogive 
with a 0.05A thick material coating containing the properties t T — 3 — 
j5,H r = 1.5 - j 0.5 
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1 x 4A Coated Conducting Ogive (E-pol) 



1 x 4A Coated Conducting Ogive (H-pol) 



Angle [deg] 


Figure 4.9: E z and H z backscatter echowidth of a 1 x 4A perfectly conducting ogive 
with a .05A thick material coating containing the properties t T = 3 — 
j 5, = 1.5 — j'0.5 


G-'j L 
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4.6 Summary 

The scattering from targets surrounded by ogival and circular boundaries has 
been presented. The finite element method produces the usual sparse sub-matrix, 
while a discrete version of the boundary integral results in a dense sub-matrix. The 
mathematical boundary enclosing the scattering structure may be judiciously cho- 
sen such that the boundary integrals are convolutional. As a result, they become 
amenable to evaluation via the FFT and leads to an O(N) storage requirement. 
Among the circular and ogival boundaries considered, the circular boundary satisfies 
the above requirements. The ogival boundary results in convolutions only when the 
source and observation points are along the same arc, while the non-convolutional 
cross- terms must be stored efficiently to guarantee the required storage requirement. 
When considering circular and ogival structures, the associated circular and ogi- 
val boundaries are usually conformal to the structure, thus providing an additional 
reduction in the number of unknowns. 

To validate the method and associated computer code, scattering patterns of 
several circular and ogival structures were given and compared with data generated 
by proven methods. 



CHAPTER V 


The Elimination of Boundary Integral Interior 
Resonances in Two-dimensional Finite Element — 
Boundary Integral Formulations 

5.1 Introduction 

The interior resonance corruption of boundary integral solutions for scattering 
computations is well known, and its treatment has been a subject of research for the 
last two decades. Methods based on the “complexification” of the wavenumber [15], 
the overspecification of the boundary conditions [16], [17], and the linear combina- 
tions of integral equations [18], [19] have been proposed, while others have focused 
on the solution technique rather than the system formulation [20]. Not surprisingly, 
when the boundary integral equation is used to terminate the finite element mesh, 
the interior resonance corruption persists and this has restricted the application of an 
otherwise promising method for large scale computations of highly inhomogeneous 
structures. 

To demonstrate the seriousness of the problem, Fig. 5.1a shows the backscatter 
echo width of a circular conducting cylinder of radius a c for TM plane wave incidence. 
The mesh is terminated on a circle of radius a Q = 1.01a c on which the boundary 
integral equation is applied. The results, displayed as a function of ka Q , were obtained 
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via the FE-BI described in chapter IV. The unknown in this implementation is 
the total field, and as seen the solution is severely corrupted, a difficulty which 
persists for other scatterers as well. (For reference, the locations of the interior 
resonant frequencies are displayed by the vertical lines at the bottom of the figure.) 
The corruption is further evidenced in the near zone scattered field plotted in Fig. 
5.1b, obtained by subtracting the incident from the total field generated via the 
aforementioned FE-BI method at ka 0 = 23.586. To render the FE-BI method 
robust at all frequencies, it is thus essential to remove the problem associated with the 
interior resonances. However, employing traditional methods such as those described 
in [2] - [6] requires either significant modifications to the original FE-BI formulation 
or substantial computing time, thus affecting the efficiency and performance of the 
method. 

Further, for the “complexification” scheme proposed in [15] to be effective, we 
verified that three different computations are required for each incident angle when 
combined with the total field FE-BI system in (4.60). That is, the total field FE-BI 
solution must be repeated three times with different complex propagation constants, 
all slightly perturbed from the free space wavenumber. Fig. 5.2, the counterpart to 
Fig. 5.1, demonstrates how the amplitude of the field quantity varies for a = 1 — 
iO.OOl and a = 1— y O.005. In the following section we present a simple modification to 
the FE-BI formulation which renders it relatively immune to the resonance problem 
without the need to repeat the solution. 

5.2 Method 

The proposed approach for eliminating the failure of the FE-BI method at in- 
terior resonant frequencies is based on the observation that the specification of the 
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tangential electric or magnetic field excitation alone at the integration boundary is 
not sufficient to yield a unique solution [19]. Therefore, we move the excitation away 
from the integration boundary by employing the scattered field as the working vari- 
able. This is accomplished by first writing the total fields everywhere in space as a 
superposition of the incident field <f> 1 and the scattered field <f> s as 

<f>= <t>' + <t> 3 (5.1) 


The boundary integral in (4.5) may be expressed entirely in terms of the scattered 
field as 


**(?)- -j£ {©(?.?.) 


dn. 


- 4>'{P« ) 


~ *‘(P . ) 


d 


dn. 


G{p,p a ) 


\dU 


(5.2) 


which differs in form from (4.5) by the incident field term. The excitation is instead 
associated with the finite element portion of the system on conduction surfaces and 
material interfaces. This becomes clear after substituting (5.1) into (4.10) to obtain 

<m e + [ w'ti’dr 

Jr* 

* 

dtfdW? 1 

+ 




d<f> 3 dW? d<f> 3 dW? 


dx dx dy dy 


[ [ox dx 


+ k 2 0 v<f>*W' | 
dpd W t e 


+ kl vpwA 


dn e 


(5.3) 


dy dy 

where the quantity ha5 been left in terms of the total field to ensure tangential 
field continuity between adjacent elements which may contain different materials. 


The expressions for TM ( <f > * = El) and TE (< p 3 = Hi) differ by the application of 
the boundary condition on perfectly conducting surfaces. For TE incidence, xj> = 0 
and the contour integral contribution on this portion of the path r e disappears. For 
TM incidence, the condition <f> 3 = —<f>' is applied after assembly, and this results in the 
elimination of the associated contour integral. Thus, in following the discretization 
procedure outlined in section 4.2.1, we obtain the system 

3 3 3 

j=i *=i j=i 


(5.4) 
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where 


u ‘ = Jj {-« 


dp dNf dp dNf' 

dx dx dy dy 


'JV'j 


+ klvpNf !> dn 


(5.5) 


Note that a e t] and b are unchanged. 

Assembling the final system and applying the appropriate boundary conditions, 
we have 


A a a A a i 0 B aa 


1 


-C/a 

Aia Aji An 0 


<t>’i 


-u, 

O 

*y 

"a 

O 




-Ud 

Maa 0 0 —L aa 




0 


for TE and 


Aaa A a l 0 — B aa 


K 


-Ua 

Aia An Aid 0 


ft 


- Ui 

0 0/0 





Maa 0 0 -Laa 


r a _ 


0 


(5.6) 


(5.7) 


for the TM case. Clearly, now, the excitation is present entirely in the FE region of 
the system, as opposed to the BI portion as was the case in the total field formulation 
of chapter IV. 

The system is then solved with the introduction of a small loss in the propagation 
constant appearing in the Green’s function in (5.2) accomplished by replacing k with 
ah, where a = 1 — jS. This substantially improved the convergence of the employed 
biconjugate gradient solver for the cases considered. However, in contrast to the 
complexification scheme employed with the total field formulation, the scattered 
field solution in (5.6) and (5.7) is relatively insensitive to 6 (provided, of course, <5 is 
very small). 



91 


5.3 Results 

To demonstrate the effectiveness, efficiency and accuracy of the proposed method, 
let us reconsider the problem of scattering by a circular cylinder via the scattered field 
FE-BI formulation. As seen in Fig. 5.3a, the far field is no longer corrupted by the 
fictitious interior resonances and the same holds for the near zone field as displayed in 
Fig. 5.3b. The results, shown for a = 1 — jO.OOl and a = 1 — j'0.005, are also seen not 
to deviate from the series solution, although the convergence rate of the solver varied 
significantly as a function of a. For example, the unperturbed (no complexification, 
i.e., q = 1) scattered field formulation converged in approximately 0.15N iterations 
over the frequency band considered in Fig. 5.3a, where N denotes the unknown 
count. For a = 1 — jO.OOl, the solution converged in 0.13N iterations whereas for 
Q = 1 — jO.005, convergence was achieved in approximately 0.08N iterations. Note 
that this fast convergence rate is due in part to the fact that the discrete boundary 
integral occupies a substantial portion of the total FE-BI system. For complex 
structures, this may not be the case and the affect is expected to be less pronounced. 

Also, we consider the TM illumination of a 5.256A x 5.256A metallic square cylin- 
der. For the implementation of the scattered field FE-BI formulation, the boundary 
integral was enforced on a circle of radius a 0 = 3.754A so that ka Q is the fifth zero 
of the Bessel function of order 6. With the incident field direction normal to one of 
rectangle’s faces, Fig. 5.4 depicts the corresponding bistatic scattered field obtained 
via the total and scattered field FE-BI formulation with a = 1 — j0.005. Clearly, 
the pattern based on the scattered field FE-BI formulation agrees everywhere with 
the moment method data. In contrast, the results based on the total field FE-BI 
formulation are substantially in error. 
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Angle [Deg] 

(b) 


Figure 5.1: Comparison of the far zone and near zone fields for TM plane wave 
incidence on a circular metallic cylinder as computed by the total field 
FE-BI method and the eigenfunction series, (a) backscatter echo width 
vs. ka 0 — the lines over the horizontal axis correspond to the eigenvalues 
of a circular conducting waveguide, (b) magnitude of the TM scattered 
field on the enclosure at the resonant frequency ka 0 = 23.586 
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(a) 



Angle [Deg] 

(b) 


Figure 5.2: Far and near zone fields for TM plane wave incidence on a circular metal- 
lic cylinder as computed by the total field FE-BI method with k replaced 
by kot in the BI equation, (a) backscatter echo width vs. ka 0 . (b) magni- 
tude of the TM scattered field on the enclosure at the resonant frequency 
ka 0 = 23.586 
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Angle [Deg] 
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Figure 5.3: Far and near zone fields for TM plane wave incidence on a circular metal- 
lic cylinder as computed by the scattered field FE-BI method with k 
replaced by ka in the BI equation, (a) backscatter echo width vs. ka 0 . 
(b) magnitude of the TM scattered field on the enclosure at the resonant 
frequency ka 0 = 23.586 





Bistatic Echo Width [dB] 
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Angle [deg] 


Figure 5.4: Bistatic echo width for TM plane wave incidence on a 5.256A x 5.256A 
metallic square cylinder with the boundary integral place at a radius 
a 0 = 3.754A 
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5.4 Summary 

A method was presented for the elimination of the interior resonance corruption of 
the boundary integral in the FE-BI formulation. By expressing the system in terms 
of the scattered field and employing a complex wavenumber in the boundary integral, 
the effect of resonances were removed. This implementation of the method was shown 
to be superior to that associated with the total field formulation presented in chapter 
IV in that only one sample per frequency was needed. Though the development 
was implemented for the two-dimensional FE-BI formulation, it is applicable to the 
three-dimensional one as well, as seen in chapter VI. 



CHAPTER VI 


A Finite Element - Boundary Integral 
Formulation for Axially Symmetric Structures 

6.1 Introduction 

A finite element - boundary integral approach is applied to the case of axially 
symmetric structures. The method follows the same procedure outlined in section 
2.2. In this implementation, the coupled potential equations [21] are discretized via 
the usual finite element method, and the resulting system is augmented by a discrete 
form of the Stratton-Chu equations [22]. The storage reduction associated with the 
boundary integral is achieved by exploiting matrix symmetries and the final system 
is computed by employing a conjugate gradient solver. 

In this chapter, the formulation for the FE-BI is described for axially symmetric 
scatterers. The results presented demonstrates the accuracy of the method along 
with showing its limitations. 

6.2 Finite Element Formulation 

In this section, we derive the analytical coupled azimuth potential (CAP) equa- 
tions [23] which are then discretized via the finite element method. A consequence 
of the formulation is a line singularity, which tends to corrupt the computed fields 
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Z 



Figure 6.1: The general axially symmetric surface with source (primed) and obser- 
vation (unprimed) points and the corresponding unit vectors 

for scattering domains incorporating lossless media. A regularization approach is 
suggested for its removal. 

6.2.1 Analytic CAP Formulation 

Maxwell’s equations in a source free region are given by [12] 


V x E(r ) = —juifiH 

(6.1) 

V x H(r) = jutE 

(6.2) 

V • D(r) = 0 

(6.3) 

V • B(r) = 0 

(6.4) 



For axially symmetric media such as that indicated in Fig. 6.1, the fields may be 
represented by a Fourier series in the cylindrical coordinate <j> as 


OO 


E(r) = J2 e rn (p,z)e jm * 

m = — co 

(6.5) 

CO 

’)«(?)= £ 
m = — oo 

(6.6) 

where 


= -L ,77(r)c-"‘rf0 

(6.7) 

(6.8) 

Up substituting (6.5) and (6.6) into (6.1) and (6.2), we obtain 


[ J 771 h mz Qg ( R^m4> )| — J^r^mp 

(6.9) 

|jrne m2 ^^(/2e m(? ! ) )| = j Hr h m p 

(6.10) 

R [ SZ ^ m P dR ^ m - r ] = J6 r (/?e m< £) 

(6.11) 

R e mp = ~j Hr{Rhm<t>) 

(6.12) 

~fi [j ^71 ^mp ( RRm<t> )| — J^r^mz 

(6.13) 

/j [j 771 Cmp Qft ( ) j = 3 Hr ^mz 

(6.14) 

with 


R = /c 0 p, Z = k 0 z 

(6.15) 

to be referred to as normalized coordinates. Substituting h mz 

of (6.14) into (6.9) 

gives 


^mp J f m TTl RCjncf)^ “f* RHt qz (-^^m^)J 

(6.16) 


where 
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and 


K 2 = ^rtr (6.18) 

Analogously, substituting e mp of (6.9) into (6.14) yields 

h-mz j fm “I" Rt T gjj(R€mii)J (6.19) 

while a similar procedure for combining (6.10) and (6.13) yields the pair 

^■mz jfm Rpr (6.20) 

^mp j fm [ m 97?(^^m^) Rt-r ~Q2 ( (6.21) 

Equations (6.16) through (6.21) may be expressed in compact form as 

4> x e m (R, Z) = jf m \m<i> x V,^ e - fi T RV t rl > h j (6.22) 

<t> x h m (R, Z ) = jf m [m<^ x V|0 A + e r i?V t V> e ] (6.23) 

$ ' &m{Ri Z) = 1pe/R (6.24) 

f> • A ro (i?, Z) = xp h /R (6.25) 

where 

= PdR + ^Wz (6.26) 


and V’e and are herein referred to as the azimuthal potentials. Rewriting (6.11) 
and (6.12) as 


RV t • {4> x h m ) = -jt r tp e (6.27) 

flV, • (^ x e m ) = (6.28) 

and then substituting (6.22) and (6.23) into them produces the CAP equations 

V ( ■ [f m (e r RV t tJ> e + m4> x V,0 k )] + = 0 

Vi • [/ m (/i r /?ViV>fc - m<£ x V<^ e )j + = o. 


f? 


(6.29) 

(6.30) 
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This system may be written in operator form as 

LV> = 0 (6.31) 

where 

Vi • [fmtrRVt] + ft mVi • \fm4> X V t ] 

L = (6.32) 

-mV, • [/ m <£ x V t ] V, • [/ m // rJ RV f ] + ^ _ 

and 

4> = tye ^] T (6.33) 

The three dimensional axially symmetric problem has, thus, been cast into a 
two-dimensional one and its discrete representation is formulated in the following 
section. 


6.2.2 Discretization of the CAP Equations 


To discretize (6.31), we first enclose the generating contours of BOR in a fictitious 
boundary T (= r a + r c + r;j) as shown in Fig. 6.2. The region interior to T is divided 
into N e triangular elements and within each element the corresponding weighted 
residual expression is 


JJ N'(R, Z) Lxl> dS' = 0 (6.34) 

S' 

where N ' is the usual shape function [5], so chosen to satisfy the Dirichlet boundary 
condition of if> e and on T*. Substituting the first of (6.31) into (6.34) gives 



(t r RV t xl> e + m4> x V t 0 fc 



dS' = 0 


and upon employing the identity 


(6.35) 


N'V t ■ A* = V, • {N'A S ) - A* ■ V t N' 


(6.36) 
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R 


Z 

Figure 6.2: Cross section of a generating surface enclosed by the fictitious boundary 
gives the expression 

JJ [v« • {N< f m (e r RV t \l> e + m<j> X 
s 

-fm (e^VtV’e + m4> x V,0*) + — y ■ V t A7 dS e = 0 (6.37) 

Furthermore, by invoking the divergence theorem, (6.37) may be written 

JJ | [-/ m (e r RV t tP c +mj>x V t ^)] - V X JV? + dS e 

+ J C ' « ' [^/m (c r RV t tJ> e + m<i> x V,^)] dl e = 0 (6.38) 

where n is the outward normal along the boundary C' of the eth element. Finally, 
on using (6.22) these may be simplified, yielding 

JJ | [~fm (trRV t tl> e + mj> x V^a)] • V ( w; + | dS e 

i S J 

- <f N'(jh mt )dl e = 0 ( 6 . 39 ) 

J O e 
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where 


with 


h m t — t • h~ 


(6.40) 


t — h x 4> 


(6.41) 


Equation (6.39) and its dual constitute a weak form of (6.31) over the eth element. 

The development thus far has employed the total potential as the working vari- 
able, but may also be expressed in terms of the scattered potential as well. To this 
end, any function V’ satisfying Maxwell’s equations may be written as a superposition 
of the incident and scattered potentials, i.e., 

& = + -0e (6.42) 

= V’j ; + Vk (6.43) 

where the superscript s denotes the scattered potential and i denotes the incident 
potential (i.e., that potential present in the absence of the scatter). Substituting 
these into (6.39), the corresponding expression in terms of the scattered potential is 

~fm (trRVttPl + x V^)] • V,A7 + dS e 

-f cc N'(jh mt )dl' = 

-fm (tr RV t il>' e + ml x V t V>;)] • V t N ,? + ds e (6.44) 

where the contour integral has been left in terms of the total field until the final 
system assembly (of equations) is performed. 

To form a discrete system of equations, the potentials in the e element are ex- 
pressed as 

r,(R,Z) = -£N’(R,Z) { 

:=i 




(6.45) 
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< K(R,z) = '£n’(r,z ) {«}■ 

J=l 


(6.46) 


and upon substituting these into (6.44) gives 


£\ fill- 

<-■ Ls' u 


N e N e ] 

RV t N'V t N' + e r -^- {#}; 


x V t N' • V t A?{^} dS ' ] 
- I N'(jh mt )dl' = {t/}' 

c 


(6.47) 


where {U}' is given below. Assuming t T and /x r are constant within the element, 
(6.47) may be written as 

1 1 £ ; Mt WYi - Wo {«};] = i W ( jh M )d r + {(/}< (6.48) 

]=l JC ' 


where 


rr T N?N e 1 

Wo- = JJ ~fm R VtN' • V t N' + dS' (6.49) 

S' L K J 

[bj'j = JJ mf m (j) x V t N' ■ V t N' dS' (6.50) 

S' 

{U}t = -JJ | [—fm (tr R + m<j> x VjV'a)] • V t A r , e + — dS e (6.51) 


(6.50) 


These constitute the equations for the eth element, and the final system is assembled 
by summing over all elements in the discrete model. 

Before pursuing the step of system assembly, explicit expressions for (6.49) and 
(6.50) may be developed by first writing them as 


w- a a^'- *//¥*• 


(6.52) 


m = -i 


(2fl e ) 2 


(6.53) 
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where the linear shape function and its coefficients are given by 


N'(R,Z) = 

of + P'Z + 7 fR 

(6.54) 

a' = 

Zj i+l n i+2 - 

(6.55) 

01 = 

^ 1+1 — ^+2 

(6.56) 

li = 

Z' + 2 - Zf +l 

(6.57) 

'■ = // 

S e 

R j ce 

k 2 R 2 — 

(6.58) 

w/ 

s e 

K 2 R 2 -m 2dS 

(6.59) 


Evaluating (6.58) and (6.59) over a triangular element yields 

'1-7EKW+M ( 6 . 60 ) 

K i-l 

= -!?(-<»)] (6.61) 
TTl »=1 

with 


l'(m) 


( kR ' -f m) log (kR' + m) 

PU1PU2 


(6.62) 


The second integral in [a]^ may be computed numerically via gaussian integration. 
The sums appearing in Iq and / 1 , however, experience problems when an element 
edge is parallel or nearly parallel to the axis of revolution (i.e., element d in Fig. 
6.3). In this case, two of the three terms have denominators which are nearly zero, 
resulting in cancellation errors upon their summation. To avoid such errors, the 
sum of two consecutive terms, i.e., J t e is carried out by first expanding the 

numerator of either term X and then performing the addition analytically. In this 
way the offending term is canceled and the final result is well behaved. 
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Figure 6.3: Typical triangular elements in the vicinity of the line singularity at Ro = 
— for real k 


Proceeding then, we note that the Taylor series expansion of xlnx about xo is 

,X - x 0 


where 


x In x = x In x 0 + (x — x 0 ) 1 + S( 


00 (_n«+i 

s (y) = L , t\ y n 

„=i n ( n + 1) 


x 0 


-) 


(6.63) 


(6.64) 


and the sum converges when |x - x 0 | < H, where H is the radius of a circle centered 
at Xo in the complex x-plane for which xlnx is analytic [24]. The series may be 
truncated in N terms when the error, given by the ratio of the Nth term to the first 
term in the series is less than some tolerance e, or more explicitly when 

2y N ~ 1 


< t 


(6.65) 


\N(N + 1)\ 

For e = 10 -7 , a fit to the nonlinear function in (6.65) for |y| 6 [0.1, 0.7] is 

N(y) = Round(e' ! 'l(9.97633el 1 'l - 5.94387)) (6.66) 

Once N is known, the series is evaluated efficiently using Horner’s rule. 

Employing the expansion (6.63) to the partial sum Ji(m) + J i+ i(m) yields 


J,(m) + I 1+1 (m) = -j 


_j_ , _ 5( .*k_)] forml 

| J_ + , + 5(3gS)] form 2 


(6.67) 
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where form 1 was obtained by expanding the numerator of T,(m) about R,+ j, while 
form 2 results by expanding the numerator of J, +1 (m) about R,. For elements well 
away from the line singularity, either form is valid. However, as kR + m — ► 0, the 
radius of convergence 71 is reduced to zero due the branch point of (/cij + m) In («;/? + 
m). The method for choosing the appropriate form may be done numerically by 
computing the value argument of S in (6.67) for each form. The one which gives 
the smallest value y € [0.1, 0.7] in S(y) is chosen for computation, since the series 
will converge the fastest for it. If y falls in the range 0.7 < y < 1, the series (6.64) 
will still converge, but very slowly as y approaches 1. If |y| > 1, the series will not 
converge and the associated form (1 or 2, or both) will be invalid. Note that neither 
of them is valid if both points lie along the line singularity (in which an infinite result 
is obtained) or if one lies above it and the other below (as in element b of Fig. 6.3). 
In the latter case, (6.60) and (6.61) must be evaluated directly, but a new mesh may 
be necessary for the former. An approach for bypassing this difficulty is discussed in 
the following section. 

Summing over all elements to obtain a solution for the entire computational 
domain our system becomes 


EE[<«« {«>;-*& {«);] = £/ N ' + e / w'( 

e=lj=l e =l- /r » e =l Jf ’ 

+ £/ Ai?yA-.)<a+E<to; (6.68) 

e=l Jlc C— 1 


Note that since Nf and h mt are individually continuous between adjacent elements 
(and the same is true for their product), the contour integrals cancel everywhere 
except on the domain boundary T. Accounting for A r f = 0 along the z-axis, this 
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system combined with its dual may be written in block matrix form as 


A‘ aa 0 A e aI -B aa 0 -B aI 




ft'). 

0 

A' c A l cI 0 —B cc -B cI 


{^e)c 


(£'}« 

A\ a A c j c A\i —Bi a —Bjc —Bu 




{£')/ 

B aa 0 B aI A£ a 0 A'j 


M}. 


(i"). 

0 

B cc B cI 0 A? c A^j 


{V’aJc 


<V'h 

B Ia B lc B n Al A 1 }, A^ 




{V'}, 


r 





T 

+ 

0 0 e„,il 

0 0 


+ 

0 0 0 -Z£,frJN?e ml 

dl 0 1 

J 

T 


(6.69) 


where 

M» = £ < Mt (6.70) 

c= 1 

ft = E A Ml, (6.71) 

e=l 

! B] k , = £ Mu (6.72) 

c= 1 

Uk = Y w: (6.73) 

e=l 

where the ij th member of the eth local element matrix is related to the kith member 
of the global matrix by a node-element connectivity transformation function p as 


k = p(e,i) (6.74) 

l = p(e,j) (6.75) 


In (6.69) the subscript / refers to group of nodes in region the fl excluding T, and the 
subscripts on A e,>J , B, U and V refer to the various regions of Q and its boundary. 
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For example, An refers to the matrix resulting from the interactions between nodes 
in region fi and on T c , whereas A cc refers to the matrix resulting from the interaction 
of the nodes on T c . Each matrix is sparse, having nonzero elements when the nodes 
share a common element. Note also that V is the dual of {/, which was defined in 


(6.51). 

Two options exist at this point for evaluating the contour integrals. They may 
either remain on the excitation side and the tangential modal fields expressed in 
terms of a condition on r o , or the tangential fields may be expressed as unknown 
scattered potentials and moved into the matrix. The second is required for the FE-BI 
formulation and is given by 

[ m a 

Ka 0 Ah -B aa 0 -Bal 0 ~C aa {r e }c {U} a 

0 / 0 0 0 0 C aa 0 {r e }i ~{4’'e}c 

A ( la A) c A)j -Bj a ~B Ic -Bn 0 0 {$[}. _ {U)i 

B aa 0 B al A£ a 0 Ah 0 0 {^}c {V} a 

0 B cc B c j 0 A& A u cI 0 0 {V h }i {V)c 

Bia Bn B u Al A% Aft 0 0 J {^} 0 [ {^}/ 

. M)c . 

,T 

+ | ES, 0 0 - E"=, fr. 0 

(6.76) 


where the boundary conditions ^6 = 0 and e mt = 0 on T c were also enforced. Addi- 
tionally, 



(6.77) 
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j= 1 

(6.78) 

«, = iRK, = EA7M,)' 

(6.79) 


J=1 


This form is now suited to augmentation by either a discrete form of a boundary 
integral equation explored in section 6.3. 

6.2.3 An Improvement to the Finite Element Formulation 

The integrand of each of [a] e tJ and [b] e {j in (6.49) and (6.50), respectively, becomes 
singular when R = ± — for real /c. Since we are only concerned about solutions 
for which m > 0 (those for m < 0 are found via symmetry considerations), only the 
positive sign is considered, or Rq = The location of the singularity is independent 
of Z and is hereafter termed w line singularity.” The line singularity intersects any 
element t containing the radius Jf?o (as seen in Fig. 6.3 for a homogeneous medium), 
or is near an element if Mixx{\RQ R* |) is small (not defined now) for the normalized 
radii of element e, R * for i £ [1,3]. The origin of this singularity has been discussed 
previously in [23] for the CAP equations. 

For elements containing i?o, [a]^ and [6]- gain an additional residue contribution, 
automatically included by the use of the logarithm of (6.62) in the sums (6.60) and 
(6.61). However, for elements oriented nearly parallel to the line singularity, as in 
element b in Fig. 6.3, the contribution to (6.60) and (6.61), and consequently [a]*, 
and [fe],j, become large. In the case of element a in Fig. 6.3, they become infinite. 
Large matrix elements are associated with the slow convergence of the conjugate 
gradient solver and inaccuracies in the resulting solution. Presented here is a way to 
avoid this problem. 
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To understand the nature of the problem, we recall (6.23) 

4> x h m (R, Z ) = jfm \m4> x 'Vii’h. + c r i?V t ^ e ] (6.80) 


the right hand side of which appears in (6.44) and ultimately in [a]^ and [6]- ; after 
discretization. Since 4> x h m must remain bounded everywhere in space (barring 
edges of perfect conductors), the right hand side of (6.80) must also. Clearly then, 
the bracketed terms involving tp e and rph must combine to cancel with the singularity 
in f m = 7 - 5 — hs — :• The finite element discretization, however, separates these 
two terms and each will individually become large when an element edge is nearly 
coincident with the line singularity. 

To regularize [a]- ; and [6]-, we first define the integral 

/ = // fm (£r RVM ■ V t N?dS e + II fm (mix V t ^) • V t N'dS e (6.81) 

s e s e 

which appears in a portion of (6.44). Each of the integrals in (6.81) must be regu- 
larized, and in doing this we must subtract and add a term from each of the form 


h(Z,Ro) 


(6.82) 


kR — m 

where h(Z,Ro) is the factor of the integrand which is well behaved at R = Ro- 
Applying this technique to (6.81) gives 


/= // 1 [ t r RVtV t . Vti ye_ £r ( Z, W Ylftl L ^ 1 .^ t N^(Z,R p) 

JJ kR — m [ kR + m 1 kRq + m 

+//; 


iR — m 


* * V,iV? - t x YMh M . V,A7(Z, Ro) 


kR + m 


kRq + m 


dS e 


dS e 


t[Z _ ', i?o)fioV t ^(Z, fip) + mcf) x V t ^(Z, .Ro) 
/cRp + m 


V t N'(Z,Ro)dS e 

(6.83) 


where /i(Z, Ro) for each integral is clear by the argument (Z, Ro)- Where not explic- 
itly shown, all functions exhibit a ( Z, R ) dependency. 
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The key to implementing this method is the elimination of the third integral in 
(6.83). Note that the numerator in the integral of this terms is simply the bracketed 
factor in (6.80) evaluated at R = Rq. As mentioned before, <j> x~h m (Ro,Z) must be 
finite, the bracketed term is zero at R = Rq. Thus, 


m<t> x V t MZ,R o) = — e r (Z, Ro)RoV t ip e (Z, Rq) 


The third integral of (6.83) vanishes when the expression is substituted into its 
integrand. 

After following the usual discretization of the azimuthal potentials, (6.83) may 
be written in discrete form as 


X/ kR - m l kR + m 


Rq)RoVjNJ{Z, Rq) 
kRq + m 


- V t 7V, 


1{Z,R«)\ 


JT X J KH — m kR + m 


<h^NJ{Z 1 R q) 

kRq + m 


■V t N, 


!(z,R o)]. 


(6.85) 


For linear shape functions, the gradients are constant (as well as e r ) and are factored 


out to yield 


7 = ■ v ‘ N ' ds ' 

+ pwjj an(^ +w) » * • v ' N ‘ ds ' 


( 6 . 86 ) 


which are clearly regular when R = f . Thus, the matrix elements in (6.49) and 
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(6.50), may now be written 


[< = //[- 

S' 

W = II 


1 N'N'l 

VN' ■ V,N' + — — L 

2 k(kR + m) 3 1 R . 


dS' 


1 


gi 2 (kR -I- m) 


<t> x v,iv; • v t N'dS' 


(6.87) 

( 6 . 88 ) 


These were implemented, but did not work very well. Note that since the potentials 
V>e(Z, Ro) and ^(Z, Ro) were expanded in terms of linear functions, their respective 
gradients are are constant everywhere in the element independent of R and conse- 
quently independent of Ro. Since the derivatives Vip e and V^v, are constant over the 
element, this approach to regularization is of low order. To improve the accuracy, the 
functional dependency of this on R must be increased and for the existing discretiza- 
tion scheme, this may be done by refining the mesh in the vicinity of the singularity. 
This, however, is not a realistic option for frequency sweeps for scatterers containing 
lossless materials. Using higher order shape functions or employing the modal field 
as the working variable (i.e., V^> e = V(i2e m ^) = KVe m $ + e m< pV R) would alleviate 
this problem. It is interesting to note that given tp e = Re m a linear variation in 
Z of e m ^, corresponds to the same in tp e . However, a linear variation of e m ^, in R 
corresponds to a quadratic variation of the same in ip e . But since t/> e is expressed in 
terms of linear shape functions here, the radial behavior is lost. 

In summary, then, we have developed a method by which to alleviate the problem 
associated the line singularity at Ro = ^ for real k. This not only eliminates the 
need for a residue contribution to [g]^ and [&]-, but also prevents them from ever 
becoming infinite. In addition, a direct consequence is that these element matrices 
are now purely real for scatterers containing lossless materials, which is not the case 
for the traditional formulation presented in section 6.2.2. 
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6.3 Boundary Integral Formulation 

In the previous section, a weak form of the wave equation was developed for 
scattering bodies in fl and bounded by T = T a + T c + T 2 . The boundary conditions 
on r c and have already been included in the system (6.76). The resulting discrete 
system remains incomplete, however, since the boundary conditions (which provide 
a relationship between the tangential E and H fields) on the exterior boundary r„ 
remain unspecified. The Stratton-Chu integral equation (S-C) provides the necessary 
relationship for source and integration points along r o , when the integral equation 
is expressed in terms of fields tangent to the surface of revolution. By employing a 
Fourier series expansion of the S-C and thereafter discretizing it on r a , the resulting 
system provides the needed equations for solving (6.76). 

As a preliminary step , the electric and magnetic fields in the unbounded region 
are represented by 

E(r) = F(r)+F(r) (6.89) 

77(r) = 77* (r) + 7T (r) (6.90) 

where E (r) and H (r) are the incident fields and the scattered fields are given by 
the Stratton-Chu equations [22] 

F(f) = ft X Wff‘(r')]s(f,f') + T^[n' ■ V' x ,„F*(r')] Vg[r,f) 

s 

+ (n x F(r')] X V' S (r,r') W (6.91) 
9o H‘(r ) = ft X £ , (f')] 9 (r, r“) - ■ V X F^)] V' S (r, ?<) 

s 

+ [n' x V7V)] x VV(r,f / )|d5' (6.92) 
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where r* and r are the source and observation points, respectively and 

p -j*o|r-r'| 

g(ry) = 


(6.93) 


47r|r — f*\ 

is the free space Green’s function. In this form, all field quantities are tangent to 
S'. Since in the development of the discrete boundary integral system the source 
and field points reside on the surface S', it is convenient to remove the singularity 
at r = r 7 by expressing the integrals in (6.91) and (6.92) in terms of their respective 
principal values as 

|F(r) = j-jfco [«' x Vo H a (r)} g( r,f) + jj- [n 7 • V' x t? 0 jr(F')] V'^r,? 7 ) 

s 

+ [n # x r^f 7 )] x V'p^rOjdS' (6.94) 
5 » 7 o H’(r) = [n 7 x ^(f 7 )] g^f 1 ) - jj- [n’ • V x Ffr 7 )] Vg(r, f 7 ) 

s 

+ [n 7 x T/o H S (r)] x Vg(r, r 7 ) W(6.95) 


Looking back to (6.76), it is clear that two additional equations relating the 
unknown potentials ip‘, x(}* t and are necessary to form a complete system. 
This may be achieved, for example, by using the i component of the modal form of 
(6.94) and the (j> component of the modal form of (6.95). In fact, many combinations 
are possible but for simplicity and symmetry, the component of (6.94) and (6.95) 
is considered below. 


6.3.1 Derivation of the Modal Boundary Integral Equation 

In a fashion analogous to the finite element formulation, the development is pro- 
vided explicitly for (6.94) and duality is employed to obtain the corresponding ex- 
pression for (6.95). First, consider the general surface of revolution indicated in Fig. 

6.1 whose tangential unit vectors are denoted by 4> and i. The unit vector t subtends 
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an angle v with the z-axis, and i' subtends an angle v' with the z-axis. With reference 
to Fig. 6.1, the various unit vectors are given by 


h = x cos v cos (f> -f y cos v sin <j> — z sin v 

(6.96) 

4> = —x sin <j> + y cos <f> 

(6.97) 

t = x sin v cos <f> + y sin v sin <j> + z cos v 

(6.98) 

x = t sin v cos <f> -f* h cos v cos <)> — <j> sin <f> 

(6.99) 

y — t sin v sin <f> + n cos v sin <f> + <j> cos <f> 

(6.100) 

z = i cos v — h sin v 

(6.101) 


Expressing the primed unit vectors in terms of the ones results in 

t — t' [sin v' sin v cos ((f> — <j>') -f cos v cos v*] 

+ h' [cost/ sin i>cos(<£ - <j>) - cosusint/] + <j>' [sin v - 4>')\ (6.102) 

h = t 1 [sin v' cos v cos(<f> — <f> r ) — sin v cos v'\ 

+ h' [cos v' cos v cos(<f> - <}>') + sin v sin i/] + <}> [cosusin(<£ - <}>')] (6.103) 

<f> = —p'sm(4> - <f>') + <j> cos(<j> - tf>) (6.104) 

which clearly reveal the dependency of the unit vectors on the azimuthal variation 
<j> — <f>'. Further, the following identities can be shown to hold: 

j>-V'g = (6.105) 

x go 7T) = sinu'sin(<£ - <f>') - g 0 H a t cos{4> - <j>') (6.106) 

ft' • (V' x tjoTT) = ^ [-^{p'vo H^) + (6.107) 

$ • j(n' x E a ) x VV] = [z'E 3 t sin(<^ - <j>') + h'E % cos (<f> - <f>') 

+ (\> E^ cost/ sin (<f> — ^')j • V'g 


(6.108) 
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Using these in (6.94) gives 

\El{f) = ^ j ^jk 0 TioH^sinv' sin{<f> -<}>') + r) 0 Ht cos{<f> - 4>')\g{r,f) 

~ j^[-MpW) + ^(VoH;)]h^9(ry)+ {z’E' t sm(<t> -<(>') 

+ n El cos(4> — <f>') + 4>'El cos v' sin(</> — <£')] • Vg^j p'dtp'dt' 

(6.109) 

and by carrying out the derivatives of the Green’s functions, we find 
\El(r) = ^ j ^ |jfc 0 7/ 0 ^sinn , sin(<^ - <f>) + rjoHj cos{<j> - <j))^g(r,r) 

- jk 0 \-Mp^ h d + sin (^ - ^R^dk 

+ f— £/(z — z') sin(<^> — <f>) + El{cos(<t> — <f>')[(z — z') sin v — (p — p) cos v'\ 

- 2pcosv' sm 2 ( - - — (6.110) 

in which 

Ro = \fp 2 + p' 2 - 2pp r cos(fli> — <£') 4- (2 — z') 2 (6.111) 

The integrand of (6.110) is expressed explicitly in terms of <f> and <j>' and is now 
suitable for harmonic decomposition. 

To generate the corresponding modal integral equations in terms of the Fourier 
coefficients, the fields and Green’s function may be expanded as 

OO 

E’(r)= £ (6.112) 

m= — 00 

OO 

-lotf'(?) = £ £.(/>■*)«*“* (6113) 

m= — 00 
00 

S li Hr,r‘) = £ j (6-114) 

n= — 00 

where 

e’ m (p,z) = f_T(p,u,z)e-‘ m -du (6.115) 
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h' m (p,z) = 
gM(p,p',z,z') = 

9 ln\p,p\ z , z ') = 

9 ( rt\p,p\z,z') = 

^‘{PiP'iZ,?') = 
9n Y (P, P',Z,Z') = 
9n ] '(P, P',Z,Z') = 


D ■ f_ v H 3 (p,u,z)e-”™du 
1 [* e - jko H 

— —cos (nu)du 

x J 0 4 nR 

1 f r e-^fi 

— 4- cosu — cos(nu)du 

X j 0 4x72 

j ri r g—jkoR 

j- sinw =r- sin(nuWu 

X J o 4x7? 

1 r 1 d 9 , v , 

— r^- -h ——=? coslnujdu 

xtfJoRdR V ' 

1 1 d# . , 

— rx f- COSU-^r—sr cos(nu)du 

xklJ 0 7?d7? v ' 

j r ■ i d 9 . , \ i 

rx 7- sm u-=r— sr sinlnuiat/ 

^o-/o 7?d7? v ' 


in which 


■ft = P 2 + p' 2 — 2pp> COS « + (2 — 2') 2 


Substituting these into (6.110) yields 


5 £ <*(/>. £ E eW rfA^i-Vf' + Of 1 ] 

n= — 00 n=-oom=-oo ^ 1 a y 0 J 

~ e mt( z - z> ) k l 9n Y + e a m<j) k 2 (p' cos v'gW 
- p cos + (2 - 2 / )sinr/^ 1) ')} e i{ ' m ~ n)4 ‘' p' d<i> dt' (6.124) 

where t is the parameter along the contour T a and increases in the clockwise direction 
as shown in Fig. 6.2. Further, upon multiplying each side by e~ Jlp<t> and integrating 
over (0, 27r) to extract the mth harmonic equation gives 

l e mt(p, z ) = 2?r ^{i^sinu^ - j^\p'h a m<t> )gW + (/*««) [rf* + i m 9m Y ] 

+ ^ m<t> k Q (p' cos v'gW - pcosv'g m + (2 - z ') sin v'g^ y ) 

~ e mt( z ~ z ')kog£ y }k 0 pdt' (6.125) 
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after also combining similar terms and where we have used the identity 



e Hm-n )4>> d( j > > 


2i r m = n 

< 

0 otherwise 


Introducing the normalized coordinates 


(6.126) 


R = hop R! = k 0 p' 
Z = koz Z' = k 0 z' 


d , d_ 

• ko dr' 


dt 


(6.127) 


in (6.125) and replacing the field quantities with the azimuthal potentials yields 


TR^e = [>«»«'*£>] - (^h) + Kt + J™g£ Y 

+ [# cos v '9 ( m - R cos vg m + (Z - Z') sin 


+ 


[j(Z - Z')gW]}dT f 


(6.128) 


This equation and its dual form a pair of integral equations to be imposed at the 
mesh termination boundary. Their discretization is considered next. 


6.3.2 Discretization of the Modal Boundary Integral Equation 

In discretizing the modal boundary integral in (6.128), the contour T a is first 
divided into N a elements each of length A e . Within the eth observation element, the 
parametric representation 


R = R{ + rsinu* (6.129) 

Z = Z\ + r cosv e (6.130) 


is adopted as shown in Fig. 6.4 and is consistent with the required counterclock- 
wise path traversal of T a . Likewise, within the source element e', the parametric 
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Figure 6.4: A typical pair of source (e^ and observation (e) elements associated with 
the discretization of the boundary integral equation 

representation 


K = R{ + r' sin v e ' (6.131) 

Z' = Zi + t 1 cos v e ' (6.132) 



are linear parameterized functions of T , . (Note that the linear basis functions are 
employed in the finite element discretization.) Employing the method of weighted 
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residuals over the eth observation element and Galerkin’s technique to (6.128) gives 

Na 2 ; 4 ' N e N e O- N ° 2 ,A* ,A«' r 

, * 2n * rA . re . re 'f. 


( 2 ) 

/m 


[;#<£''] + <«,};' [»s» + w«'] 

+ {«)i AtfJVflfl' «*»'«£>■ - Rcosv'g' m + (Z - Z')sin ..'si, 11 '] 

+ {«}>' W?<[j(Z - Z')jW] }*•'*• = 0 (6.136) 

In deriving the first term in equation (6.136), it is useful to note that 4>Z(t) may be 
written 

r,{r)= L r,(r')S(r - r')ir' (6.137) 

r a 

where 6 is the Dirac delta function. Substituting (6.133) into (6.137) yields 

N 2 

«!■) - E EWtf L N’’(t')6{t - r')dr' 

e'=l ; = 1 • /r « 

= E EW,');'<M 

e'=l j=l 

Taking the inner product of (6.138) with yields 

V^^r./aie' / (r) 


(6.138) 


/ r fw=ED«; 7 r 

• /I ° n e'=lj=l ' /r 




dr 


(6.139) 


as expected. 

The system (6.136) may be written in compact matrix notation 

N„ 2 

EE 

e' = l j - 1 1 

where we easily see that 


E E Mfttnj' + + N“' + [m, 


= 0(6.140) 


: cos i>'^ 


, r&‘ N'N e ‘ ‘Itt r&* /-A*' 

MS—il /„ *r*r[*w*E ), -Ji. 

+(Z — Z') sin v'g^ j dr' dr ( 6 . 141 ) 
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!«,]*■ = [;«>'] }<*t'<;t(6.142) 

[i.]”' = Y o jf'/f N;Nf[j(Z - Z‘)g^}dr'dr (6.143) 

[M]"’ = / 0 NtNf\ g V+jm£r\dT'dT (6.144) 


(6.144) 


Each matrix of the form [>4]" can be termed “element interaction matrix” (after 
the finite element term “element matrix” for the case e = e'), since it represents 
the interaction between the source element e' and the observation element e for 


iij € 1,2 of the local weighting and expansion functions, respective!}'. The numerical 


evaluation of the integrals in (6.141) - (6.144) is performed by breaking the integral 
into two integrals with continuous integrands and performing midpoint integration 
for each of them. This applies to the integration in r' as well as for r. When the 


source and observation points coincide, an average value is computed by moving the 
observation point away from the sub-cell center. This has been shown to work 
well in [25], and avoids the need for elaborate self-cell evaluation techniques. In all 
cases, Gaussian quadrature is employed for the <f>' integration in (6.117) - (6.122). 

The local “interaction element matrices’ may be assembled to form a global 
boundary integral system in a fashion analogous to the finite element method by 
first summing (6.140) over all observation elements e as 



[l aM ml + mvimf + 


= 0 
(6.145) 


where the local element subscripts i,j have been replaced by their global counterparts 
k, l and are related by the node-element connectivity transformation function q as 


k = q{e,i ) 

l = <l(e',j) 


(6.146) 

(6.147) 
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A(k,l) = Ok,l€ {1,7V 0 } 

Do e = 1 to N a 

Do e' = 1 to N a 
Do i= 1 to 2 

Do j — 1 to 2 
* = q(e,i) 

l = 9(e',i) 

Form [A]£ 

A(M) = + M"' 

End Do 
End Do 
End Do 
End Do 

Table 6.1: An algorithm for the boundary integral system assembly for a generic 
element interaction matrix [A]*J 

In Table 6.1 is illustrated an algorithm where a generic matrix [A^ (where A 
represents any of Mj,, L t or M t ) is assembled to form a global system. This 
algorithm differs from that of the finite element method, in that the loop in e' is 
eliminated. 

The final FE-BI system is formed by augmenting finite element system with 
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(6.145) and its dual giving 


A c 

10 

0 

Vi 

-B aa 

0 

— Bal 

0 

-C aa 


' «}. " 


{U}a 

0 

I 

0 

0 

0 

0 

Caa 

0 


{«. 


-ii’Dc 

Va 

Vc 

Vi 

— Bla 

—Bic 

-B n 

0 

0 


{ </’,•); 


{U}i 

B aa 

0 

Bal 

Ka 

0 

Ki 

0 

0 


m. 


{V}a 

0 

B cc 

B cl 

0 

Vc 

Ki 

0 

0 


{«}< 


{V}c 

Bla 

Bic 

Bn 

a 

VC 

Vi 

0 

0 


{«}/ 


{V)i 

L<t, 

M* 

0 

0 

0 

0 

L t 

M t 


{«.}. 


0 

_ -M* 

L* 

0 

0 

0 

0 

—M t 

L t 


r 

o 

1 


0 

- 


+ 


Efi, Jr. WhLMI 0 0 -&Jr.N;Ue'„,)dl 0 0 


i T 


(6.148) 


which can be 


solved by the conjugate gradient method. Before considering different 


solutions of this system, we must develop the harmonic coefficients of various sources 


of excitation. This is considered next. 


6.4 Sources of Electromagnetic Radiation 

The modal forms for a plane wave excitation and an electric dipole source are 
developed in this section. The former is employed for scattering problems while 
the latter may be used for both scattering and radiation. For both cases being 
considered, the resulting expressions for the azimuthal potentials are used in the 
expressions developed in section 6.2. 
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6.4.1 A Plane Wave Excitation 

A plane wave incident at angles (O', cf>') and observed at r = (r, <f>, z) (see Fig. 6.1 
) has the form 

= Pe-* 0 ' 7 ( 6 . 149 ) 

u e (9\ p, <t>, z) = 0'e~^ 7 ( 6 . 150 ) 

A , * . 

where <f>' is perpendicular to the plane of incidence and O' is in the plane of incidence. 
Expressions (6.149) and (6.150) may be written as explicit functions of <j> by first 
writing the argument of the exponential function as 

k 0 -f = k 0 r(-z ■ f) 

= —ko [psin O' cos(<f> — <f>') + zcos O' j (6.151) 

since 

f — x sin 0cos <f> + y sin 6 sin <f> + z cos 0 (6.152) 

z = x sin O' cos <f>' + j/ sin O' sin <f>' + z cos O' (6.153) 

Also, making use of the unit vector transformations 


4>' = — x sin O' + y cos O' 

(6.154) 

O' = x cos O' cos cf>' + y cos O' sin <f>' — z sin O' 

(6.155) 

x = p cos 4> — <j> sin <j> 

(6.156) 

y = p sin <f> + <f> cos <f> 

(6.157) 

z = z 

(6.158) 


(6.149) and (6.150) can be rewritten as 

u^(0'\ p,4> — <j>', z) = [psin(<^ — <f>') + 4>cos{<j>- e jko[ps '•«* <***-**)+*«**•] 
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(6.159) 


Ug(0'; P,4> — ft, z) = p cos 0' cos{4 > — ft) — j>cos0' s\n(<p — ft) - z sin 0' 


gjko[p&in$ % cos(<t>—<t>')+z cos0‘] 

These may be expanded into a Fourier series in the 
(6.149) and (6.150) as 


(6.160) 


parameter (<f> — ft) by first writing 


OO 

u^ i \P^-<f>\z)= £ u m 4{0'-p,zy m(<t> -^ (6.161) 

m=— oo 

u e (0'-,p,4>- ft,z) = u me (6» , '; / 9,2)e jm ^- 0,) (6.162) 

m=— oo 

Each component of (6.159) and (6.160) may be expressed in terms of one of the 
functions 


/(0« ; p , (f, _ «£') = e i*o/»«nf cos (*-*•) 


(6.163) 


fc{6';p, 4> - <t>') = cos(</> - ft)f{0 K , p,<i>- ft) (6.164) 

- <t> x ) = sin (<j>— ft)f(0';p,<f>- <j>') (6.165) 


Expanding each of these into a Fourier series in ( <j) — ft) and using the even/odd 
properties 


S {4>) = s(-<f>) <=> s m (u) = s_ m (u) for f,f c 
s{<j>) = -s(-<f>) <=* s m (u) = s_ m (u) for f 3 

of the functions / in (6.163) - (6.165) may be expanded as 


f(0'-,p,<f>- ft) = f o (0\p) + 2 £ / m (0>)cos[m(^ - ft)) 

m=l 

oo 

/c(0‘; PA- P) = fco(0\ />) + 2 £ fcm(0\ p) cos[m(^> - ft)) 

m= 1 


Js{.9 x \P,4> - ft) = 2; ^ f am (0',p)s\n[m(<j>- ft)] 

m= 1 


(6.166) 


(6.167) 

(6.168) 


(6.169) 
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with the corresponding Fourier coefficients 


fm 

II 

1 ft 

1 k 0 p sin 6' cos u \J„ 

— 1 e J cosimujdu 

(6.170) 


it Jo 

f cm 

II 

- r cosue jk0p *‘ ne ' cosu cos(mu)du 

(6.171) 


it Jo 

f am 


-- r sin ue ]kop * ine ' cosu sin(mu)du 
7 r Jo 

(6.172) 


Comparing these to the Bessel function identities 

= - r e-’ 0cosx cos(mx)dx (6.173) 

7 r Jo 

— / cos ic ,0COSI cos(mxWx (6.174) 

7 T Jo 

Tfl 7 f ^ 

- — j m Jm{P) = / sin xe^ COST s\n(mx)dx (6.175) 

P Tt Jo 

where the last two are derived from the first by differentiation with respect to (3 and 
integration by parts respectively, we conclude that 


fm{0',p) = j m J m {k 0 p sin O') 

(6.176) 

fcm{9\p) =j m - 1 J’ m (k 0 p sin 0<) 

(6.177) 

/«.(*>) = -t m . t U«\ P ) 
kop sin O' 

(6.178) 

Using these, along with (6.163) - (6.165) in (6.159) and (6.160), the Fourier coeffi- 
cients of (6.161) and (6.162) are found to be 

= e jk °’™ 9i [pf am (0\p) + J>U(e\p)\ 

(6.179) 


u me (9';p,z) = e i k ° zcos9 ‘ [pf^O'.p) cos O' - 4>f sm {0\p) cos O' - z f m (0\ p) sin O' 


or, using (6.166), 

OO 

f>m{0',p) sin[m(^> — <f>')\ 

m=l 


u<t>{9'\p, <f> — <j>' , z) = e- , * 0 * CO8 *' 


(6.180) 
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m=l 


+ <t>fco{0',p) + 4>2 ^2 fcm{0\p) cos[m(<£ - (f> l )]\ (6.181) 

u e {^\P,4>-4>\z) = e jk ° zcote ' 

l^cos^ U($\p) + 2 £ fcm{B\p)c os[m(^-^)]| 

. r °° 

4> cos 9 X 2j ^2 fsm{0',p)srn[m((f> - <£')] (6.182) 

m=l 

. r °° 1 ■) 

— z sin O' f o (0',p) + 2 £ f m (0' , p) cos[m(<t> - [ 

771 = 1 _ J 

The modal coefficients u me and u m<t> can now be used for determining the harmonic 
coefficients of the incident and magnetic fields. For TE Z incidence, we have 


e rn = (6.183) 

h'm = «mfl (6.184) 

and for T M z incidence 

e' m = (6.185) 

h m = u m<t> (6.186) 

More explicitly and in component form 

<* = /-(*>) (6.187) 

= f.m(6‘,p)cose' (6.188) 

e mi = 7>)sini) eJ**- 1 "*' (6.189) 

C = [/m(^'./>)sm 6* Cost) — /cm (#’,/>) cos 0* sin uj (6.190) 

for TE Z polarization and 


4* = ~f>m(0\p) cos 6' e jkozcoa8 ‘ 

hint = /«-(*>) e 7 * 0 *' 08 *' 


(6.191) 

(6.192) 
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<4< = [/cm /j) COS sin t; - / m (^', />) sin cos vje^ 0 ^ 089 ' (6.193) 

h'mt = /«,(*>) sin t; e jko * c °° 9 ‘ (6.194) 

for T M z polarization, where v is the angle between the z-axis and the vector t in 
the p — z plane. The corresponding potential forms are obtained by multiplying the 
previous expressions by R = k 0 p. 


6.4.2 An Electric Dipole Source 


An i-directed electric dipole of moment xAnt/k 3 located at a point 2 0 on the 
z-axis is given by [26] 


E r = 


E e = 


e \kr(kr — kzocosd) ^3 3 j 

“rT 

p-jR o 




(± + 2L-A-1 
^ + ilo > R 


[ krkzo sin 2 6 / 3 3 j 


R o 


Ro 
1 


+ 1 


( fl 2 + K - 1 )-(^ + ^- 1 ) cos9 


Ri K Ri R o 


Es = 


e~ in « r j 1 


7T + -Ao - 1 sin <^ 


Ro IRo R% 


vHr = j ■ 


o J Ro 


vHe = j 

V H* = j 
where 


.e 


R o 
-jft o 


’ j 1 1 

— + -^2 J kz 0 sin 6 sin <f> 
j 1 1 

-p. — h -pzz | (A:r — kz 0 cos 0) sin <t> 
L/to Kq} 


Ro [Ro ' Rl 

j i 

~^~[Ro + M 


( kr cos 0 — kz 0 ) cos <f> 


sin 6 cos d>(6.195) 
cos <f> (6.196) 

(6.197) 

(6.198) 

(6.199) 

( 6 . 200 ) 


Ro = k\J r 2 sin 2 0 -f (r cos 0 — z 0 ) 2 


( 6 . 201 ) 


When using this source in connection with the current implementation, the harmonic 
coefficients of the fields must first be determined. This may be done rather trivially 
by noting the identities 


sin <j) — 
cos <f> = 


2 j 

+ e -j* 


( 6 . 202 ) 


2 


(6.203) 
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Introducing these into (6.195) - (6.200) and comparing with (6.161) and (6.162) we 
find (on setting <f > * = 0) that 


D -jRo r 


&mr — 


^m.6 — 


kr(kr — kz 0 cos 9) ( 3 i 3,7 1 ^ j 1 


2/?o L 

[■ krkz 0 sin 2 0/3 

~i% T r o 


£m<t> 

kmr 


2Ro [ 


Bo R* 


(^ + %-l)-(jr + -^-i)co ,0 


e -jR o j- j 


1 


+ Tv> — 1 


2/2o [Ro 


0-jR O 


+ -irl kz 0 sin 9 


2R* Lflo J 


sin0 (6.204) 

(6.205) 

(6.206) 
(6.207) 


e j fi ° j 11 

hmB = + ~ kz ° cosd ) 


hm4> — J 


o-jR O 


1 


2/2o l/2o + 


(fcr cos 0 — kz 0 ) 


(6.208) 

(6.209) 


where m = ±1. For all other m, the coefficients are zero. Additional field components 
needed in the finite element portion of the system are e mt and h mt and are given by 


e m t = e mr cos(0 — v) — e m6 sin(0 — v ) (6.210) 

hmt = h mT cos(0 — v) — h m e sin(0 — v ) (6.211) 

where the identities 

r = t cos(9 — v) + n sin(9 — v) (6.212) 

9 = -t sin(0 - v) + t cos (9 - v ) (6.213) 

have been involved in deriving (6.210) and (6.211). 

6.5 Scattered Field Computation 


Once the modal coefficients are determined from the solution of (6.76), the goal is 
to proceed with the determination of the scattered or radiated field. In the following, 
the evaluation of these from a knowledge of the modal scattered field coefficients is 
discussed. 
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The scattered fields in the far zone are given by 

E\r) = E^{r)4> + T]oHt(r)0 
VoH’ir) = r?o H4r)4> - E4r)0 


(6.214) 

(6.215) 


These involve the <f> components of the electric and magnetic fields and these 
related to the modal field components via the Fourier series as 


are 


M 


El(r,0,(t)) = 2 j J2 e m*(M) sin(m^) 

m—l 

M 

= ho(r, 9 ) 2 53 A^(r, 9) cos(m<f>) 

m=l 

and that for the ^-polarized receiver, and by the series 

M 

&■> 4 ) = e o( r ? 0) + 2 5Z e^(r, 9) cos (m<^) 

m=l 

M 

T 1oHl(r,9,<f>) = 2j J]] hm<t>( r i Q) sin(m^) 


(6.216) 

(6.217) 


m=l 


(6.218) 

(6.219) 


for the ^-polarized receiver. For a unity amplitude incident field, these imply the 
radar cross section formulae (2.33) 


cr = lim — 

r — oo 471- 


cr<t>e = lim — 

r — oo 47T 


<?$$ = lim — 
r — 00 47T 


(Tee = lim — 

r-.oo 4^- 


M 


+ 2 5Z ^ 0 (r, 0) cos(m<t>) 

m= 1 
Af 

^00 + 2 13 ^m<*i( r i 0) cos(m^) 

m = l 
M 

2 ]£ ^m^( r > 0) sin(m<^) 

m = l 
A/ 

2 H Q)sin(m<f>) 


m= 1 


( 6 . 220 ) 

( 6 . 221 ) 

( 6 . 222 ) 

(6.223) 


where 


c ‘^ = 4 xr e* 


m0 


= 4?rr 


(6.224) 

(6.225) 
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The evaluation of e 3 m<t> and in the far zone using the near field values. 

The Stratton-Chu integral equation used for calculating the scattered fields was 
discretized in section 6.3.2. For the potential distribution computed by solving the 
system (6.76), and these values are known on any path T / passing through the 
solution region fl and which encloses the scatterer. By subdividing T/ into Nj 
contour (boundary) elements, the scattered field may be expressed in a form similar 


to (6.136) as 


where 


and in which 



n, 2 

p* = E Efntf 

e'=l j=l 

q * = e Ew*};' 

e'=l j=l 
N, 2 

ft = EY.W 

e'=l j=l 

N, 2 

Qi='£ EW.tf 

e'=l j= 1 


(6.226) 


(6.227) 

(6.228) 

(6.229) 

(6.230) 


{P*)1 

ml 


^ jf Nf [# cos v'g£ y - R cos v'g' m + (Z - Z') sin vg^' dr' 

Jo W [ J sin U ^ 2) ] “ } dT> 

%J*\f[ j{ Z-Z') 9 Zy} d r> 


{ p, 

i Qt }]' = -j^f 0 N i [5m ) + dr' 


(6.231) 

(6.232) 

(6.233) 


(6.234) 
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which must be specialized for far zone (r — ► oo) computations. In deriving expres- 
sions for the far fields, the argument of the Green’s function is replaced by its usual 
far zone approximation 

R = \J R 2 + R n — 2 RR! cos u -f- (Z — Z') 2 ~ kor — Z ' cos 0 — R f cos u sin 0 (6.235) 
where 


Z = k Q r cos 6 
R = k^r sin 0 

The Green’s function and its derivative may then be written as 


p-jko r 

n cos 6+ Ft' cos u sin $) 

9 Airr 
Idg j 

RdR~ k 0 r 9 


and when these are used in (6.118) - (6.122), we obtain 

,-jkor 


9 {1) 

jm 


4xr 


-fcm{0,R!) e 


j Z' cos 6 


*-jkor 

s "~-b^ ue ' R)e,z '"' e 


where 


f m {0,R!) = j m J m {R‘ smO) 
U{e,R') = j m - l J' m {R' smQ) 


m 


/ 4m (0,/O- R , sm6 f^^) 


(6.236) 

(6.237) 


(6.238) 

(6.239) 


(6.240) 

(6.241) 

(6.242) 

(6.243) 

(6.244) 


(6.245) 

(6.246) 

(6.247) 
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Figure 6.5: The configuration of a typical boundary element e' passing through tri- 
angular element e 


Furthermore, substituting (6.236) and (6.237) into (6.231) - (6.234) yields 


, e~* k<>T r*c' , 

i P 4>Yj ~ j 0 N i ( T, )[ cos 0 sin r'/cm — sin 0 cos v'f m ] e jZ ‘ cose dr' (6.248) 

, e -jfc ° r r A ,' 

J^T J Q NJ {t') sin v'f am e jZ ' cos 6 dr’ (6.249) 

{ Pt }) ~ -j * 2 -j— l N i Y') cos 6f am e iZ ' cose dT' (6.250) 

, p-j^or M.i 

(«<)* ~ l jr—iT' (6.251) 


^ e mi a °d h m t are not known on Ty, then they must be computed by numeri- 
cally differentiating the potentials for each triangular element traversed by T f. For 
convenience, it is assumed that the e'th boundary element passes through the eth 
triangular as indicated in Fig. 6.5 (i.e., the boundary element must begin at one 
edge of the triangle and end on another). In this case, e a mt is found by first taking 
the dot product of n' (the outward normal of T/) with (6.22) to obtain 


J e mt = fm \m t • ~ Rfl- 


(6.252) 
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Since in the eth triangular element 


we have 


where 


j= 1 


v.mj* = Ev A 7<^'); 

J = 1 


» 20' 

V,A" = '2i±E£ 
1 20 ' 


Further, on noting the identities 


t' • z' = cos v' 


h' ■ z' = — sin v' 


t' • p' = sin v' 


*/ / 
n • p = cos v 


the derivatives can be written more explicitly as 


j=i 


(6.253) 


(6.254) 


(6.255) 

(6.256) 


(6.257) 


(6.258) 

(6.259) 


A similar procedure is used for h' • V t {^} e a nd i' * an< ^ substituting these 

into (6.252) gives 


wU&MU = 




J+5 




e' l2 

'j+2 


’^-— 2 [ mi' ■ V,r, - ■ V«] (6.260) 


= -p }'7jr m , K • v,« + ■ v«j <6.26i) 


where is introduced to express the left hand side as an azimuthal potential. 

The second of these is derived by a similar procedure, or via duality. For the special 
case of Tj coinciding with T a , then the e‘ mt and h* mt are computed during the system 
solution. Then the need to compute the t— components is unnecessary. 
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6.6 Code Validation and Performance 

In this section, the numerical implementation of the FE-BI method is evaluated 
in terms of storage and computational efficiency and accuracy. Also, the method for 
eliminating the boundary integral resonances presented in chapter V is tested here 
for spherical enclosures. Finally, the scattering from several structures is presented 
for validation purposes. 

6.6.1 Storage Efficiency 

This implementation of the FE-BI method discussed in this chapter differs from 
the two-dimensional FE-BI formulation discussed earlier in that arbitrarily shaped 
mesh terminations are allowed. As a result, an increased storage efficiency is not 
obtained from the convolutional properties of the boundary integral operator, but is 
instead realized by exploiting matrix symmetries. It is, nevertheless, understood that 
the controlling factor in the storage of large FE-BI systems is the dense boundary 
integral subsystem, which grows as 0(n 2 ) for n unknowns on the boundary, while 
the interior FE system grows as O(N) for N unknowns in the interior region. 

It is possible, however, to choose an enclosure on which some of the integrals 
are convolutional. The Green’s function in (6.93) appearing in the integrand of the 
integral equations (6.91) is a function of the distance between the source and field 
points, which when expressed in cylindrical coordinates takes the form 

|r — r'l = yf? + p'2 - 2p/>' cos(^ - <£') + (z - *') 2 (6.262) 

After employing the Fourier series expansion of the field quantities and Green’s func- 
tions, the azimuthal dependency is removed and the result is given by (6.128). This 
expression is a convolution only on contours parallel to the z-axis. Consequently, 



137 


a suitable enclosure is rectangular in shape, generating a right circular cylinder in 
three-dimensional space. Though analogous to the two-dimensional rectangular en- 
closure explored in chapter III, the integrals are in the form of convolutions only 
when source and observation points lie on the portion of r a which is parallel to 
the z-axis. All remaining terms are considered “cross terms” and must be stored 
efficiently, relying on symmetry to achieve this goal. 

The discrete version of a typical boundary integral equation block matrix (of L t , 
say) has the form 

an a 12 G13 

a 2l a 22 a 23 (6.263) 

<*31 G32 033 

where the subscripts of a p9 denote the observation on side p and source on side q of 
the rectangular boundary (the z-axis is not applicable to the boundary integral). As 
a consequence of the convolution, <122 is Toeplitz in form and, effectively, only the 
first row need be stored when an FFT is employed to evaluate the associated matrix- 
vector products. For 022 to dominate the storage requirement of the BI system, the 
enclosure must be long with respect to its radius. For scatterers not satisfying this 
requirement, the reduction of <122 is not substantial. Consequently, the method was 
abandoned in favor of a general conformal boundary. 

To determine the storage requirement for a general boundary, it is first noted that 
each matrix in the BI subsystem of (6.148) need only be stored once. Additionally, 
after examining (6.141) - (6.144) M t and L t are the only symmetric matrices. Thus, 
for n nodes on the boundary the storage for the each of the symmetric matrices 
is n(n + l)/2 and that for the unsymmetric matrices is n 2 . The overall storage 
requirement of the BI matrix is then n(3n — 1). Though this continues to be 0(n 2 ), 
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the constant in front is reduced from 8 to approximately 3. For the particular case 
of a purely smooth metallic scatterer, the full storage requirement of the standard 
surface integral equation (SIE) for BOR structures is 4 n 2 , slightly less than that of 
the FE-BI system. However, for structures such as a corrugated cylinder, n becomes 
very large for the SIE but in the case of the FE-BI method, the termination boundary 
may be chosen to occupy a minimum path and thereby minimizing the associated 
storage requirement. 

In passing we note that it is possible to optimize the storage by choosing the 
termination boundary far enough from the scatterer so as to reduce the sampling the 
requirement of the boundary integral (which grows as 0(n 2 )), but without substan- 
tially increasing the storage associated with the additional elements required in the 
finite element region (which grows as O(N)). This is particularly useful for structures 
containing very thin high contrast material layers, in which case a strictly conformal 
boundary would result in a large number of unknowns on the boundary. 

For a given amount of data storage space, S [bytes], the maximum allowable 
length / of the termination boundary depends on the uniform sampling rate A as 

l = nA (6.264) 


The single precision (8 bytes) storage requirement for the discrete boundary integral 
thus becomes 


S = 24 



and then solving for the length, 


/ = A 



(6.265) 


(6.266) 


Thus, if S = 10 megabytes and A = A/20 


/ ~ 32A 


(6.267) 
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This represents an upper limit on / given the memory and sampling rate without 
regard to the finite element storage nor the overhead associated with the geometry 
and solution storage, which mutually grow as O(N). Consequently in practice, n = ^ 
must be reduced until the storage of the entire model is accomodated, a requirement 
which depends on the scattering body. 

6.6.2 Computational Efficiency and Accuracy 

The FE-BI system given by (6.148) is solved via the conjugate gradient method. 
It is well known that the rate of convergence is proportional to k 2 , where k is the 
condition number of the matrix. Clearly, those factors which influence k must be 
examined, since in part, cpu time of the solver and subsequently the accuracy of the 
solution depends on this. 

Since the FE-BI system is comprised of two incomplete subsystems, it is expected 
that the weakness from each portion contributes to a reduction in k. Among those 
due to the finite element subsystem, are the shape and size of the element, and 
also the lines of singularity at Rq = ^ for real k, discussed previously. Elongated 
elements yield larger matrix values than an equilateral triangle of the same area, and 
the line singularities produce matrix values roughly two orders of magnitude larger 
than the average matrix value. The combination of large and small matrix elements 
can undermine the condition of the matrix. 

Those factors influencing the BI portion are the element size and the total number 
of boundary elements. Elements which are very small (< 0.01 A) result in subsequent 
rows which are nearly equal and give the appearance of a singular matrix. Large 
dense systems are known to have a smaller k simply because of the number of matrix 


elements. 
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(a) 



(b) 

Figure 6.6: The mesh used for the frequency sweep of a conducting sphere is shown in 
(a) and the expanded region shows the values of ka at which the line sin- 
gu arity intersects the nodes. The frequency sweep in (b) demonstrates 
the inaccuracies associated with the line singularities 
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ka 

Figure 6.7: A sphere coated with /i r = e r = 1 — j5 is swept in frequency 

The issue of accuracy, which is closely tied to the issue of system condition, 
is examined as a function of frequency for two spheres. Each is comprised of a 
conductor of radius a with a mesh termination at a 0 = 1.02a. The second sphere 
contains dielectric material with t r = fi r = l — j5 in the finite element region. As 
seen in Fig. 6.6, there are locations where the solution has a large error with respect 
to the exact result [27], plotted as a function of the outer radius normalized to the 
free space wavenumber A:©. The vertical lines denote the regions of greatest change, 
and from (a), it is clear that the edge of an element is nearly parallel to this line. This 
problem, discussed previously in section 6.2.3, is due to the fact the matrix elements 
are becoming large and the sensitive cancellation required is not occurring accurately. 
When a lossy coating is introduced, the line-singularity is no longer present, and as 
shown in Fig. 6.7, the results follow quite well with the moment method data. 
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6.6.3 Elimination of the Internal Resonance Corruption of the Boundary 
Integral 

The corruption of the solution due to interior resonances of the boundary integral 
perseveres in the three-dimensional case. The approach taken in chapter 3 is also 
employed here. It is again shown that for the coated sphere test structure, the 
“complexification” of the wavenumber in conjunction with the scattered field FE-BI 
formulation successfully removes the resonances. Since the theory surrounding the 
method was previously presented, it will not be repeated here. 

Consider the conducting sphere of radius a enclosed in a boundary of radius 
a 0 - 1.02a and e r = 2 - j'4. Shown in Fig. 6.8 is a as a function of ka 0 of the 
outer boundary for axial incidence and observation (thus requiring the solution for 
the m = 1 mode). Clearly, the expected resonances indicated by the vertical lines 
are present and are evident in the solution for a = 1 by the spiked behavior, except 
at k 0 a o ~ 10.6 at which the resonance is shifted due to numerical errors. These are 
eliminated by setting k = ako and the results for a = 0.005, o = .007 and a = .01 
are presented in the same figure. As seen by these it is clear that the solution is 
relatively insensitive to a. 

6.6.4 Scattering from Various Test Bodies 

To validate the code, bistatic scattering patterns are considered in the four various 
sample targets. In each case the radar cross section in (2.33) is computed and the 
results are compared to a moment method solution [25, 19]. Most of the cases 
considered involve axial incidence (along the z-axis) and in this case TE Z and TM Z 
no longer have meaning. Planes defined by the electric field vector and the z-axis 
form the E-plane, and the magnetic field vector and the z-axis form the H-plane. 
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Figure 6.8: The axial backscatter cross section is displayed as a function of normal- 
ized radius k 0 a o . The structure is a conducting sphere of radius a coated 
with a dielectric with e r = 2 - j4 has an outer radius of a 0 = 1.02. 
Employing a complex wavenumber k = &oa, the resonance behavior is 
removed for the complex values in comparison to the a — 1 case 
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The patterns shown are computed in these planes. 

First we consider the scattering from a coated conducting sphere of radius 0.098m 
with a 0.02m coating with t r = 2 — j4. Fig. 6.9a demonstrates the convergence of the 
solution as more Fourier harmonics are summed together, and shown are the partial 
sums up to m=10. Fig. 6.9b shows the final converged result of the (a) (plotted on a 
scale shifted by -60 degrees to coincide with the pattern due to axial incidence) along 
with the axial incident result (requiring only m=l). Only the E-plane patterns were 
computed. A comparison to a method of moments result demonstrates agreement. 

Next we consider the scattering from a perfectly conducting sphere of radius 1A 
as indicated in Fig. 6.10. Both E-plane and H-plane patterns are included and the 
comparison with moment method results is favorable. 

In Fig. 6.11 is shown the patterns for a 2A x 0.176A perfectly conducting ogive 
with a dielectric coating of thickness 0.05A with e r = 2 — j>2, shown for an axially 
incident plane wave. Both E-plane ( TE Z ) and H-plane ( TM Z ) are shown along with 
data generated via the moment method. 

Next we consider the case of a coated sphere-capped cone frustum, the mesh of 
which is indicated in Fig. 6.12. The structure is 4A in length, the spheres are of 
radii 1A and ^ and the coating thickness ^j. Additionally, the material is dielectric 
with t r = 4 — j5. The bistatic cross section computed for this structure for an 
axially incident plane wave is indicated for T E z (or E-plane) and T M z (or H-Plane) 
incidence in Fig. 6.13. Reference data provided indicates a good agreement. 
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Angle [deg] 

(a) 



(b) 

Figure 6.9: The bistatic scattering pattern for a coated sphere, (a) Shows the sum- 
mation of modes converging to the correct solution for an incidence of 
60 degrees. The converged solution for modes 0-10 at 60 degrees (plot- 
ted on a scale shifted by -60 degrees to coincide with the pattern due to 
axial incidence) shown with the axial incident result and a comparison 
to moment method results 
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Figure 6.10: The bistatic scattering pattern for a perfectly conducting sphere of ra- 
dius 1A. Both E-plane ( TE ) and H-plane ( TM ) are shown along with 
data generated via the moment method 
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Figure 6.11: The bistatic scattering pattern for a 2A x 0.176A perfectly conducting 
with a .0.05A coating with e r = 2 — j2 for an axially incident plane 
wave. Both E-plane ( TE ) and H-plane ( TM ) are shown along with 
data generated via the moment method 



Figure 6.12: The mesh of the sphere-capped cone frustrum 
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(a) 



(b) 


Figure 6.13: The bistatic scattering pattern is shown for a plane wave axially incident 
on the coated sphere-capped cone frustum. Part (a) shows the E-plane 
pattern, while (b) demonstrates the H-plane pattern 
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6.7 Conclusions 

We have presented the finite element - boundary element formulation for axially 
symmetric bodies. The storage requirement associated with the boundary integral 
was reduced by symmetry considerations. Several patterns were generated to demon- 
strate the validity of the method, and the model shortcomings were presented as well. 
Also, a method for eliminating the problem of the line singularity inherent in the 
CAP finite element formulation was given and was shown to give real matrices for 
lossless structures. It should be noted that frequency sweeps and single mode radia- 
tion problems make the most efficient use of the employed CG solver. The method 
is not well suited, however, to problems incorporating multiple excitations or those 
requiring many Fourier modes. 



CHAPTER VII 


Conclusion 


7.1 Summary 

In this thesis, the FE-BI technique (presented in general terms in chapter II) 
was developed for two-dimensional and axially symmetric structures. The two- 
dimensional case was based on the moment-method version developed by Jin [9]. 
An FE-BI formulation for rectangular enclosures was developed in Chapter III and 
lead to simple boundary integrals some of which had convolutional form which was 
exploited for reducing the memory requirement. In chapter IV, circular and ogival 
enclosures were considered. For the circular boundary, the boundary integral was 
entirely convolutional in form and an 0(N) storage requirement was thus achieved 
at all times. It was shown that circular enclosures are consequently preferred for 
storage efficiency if the structure’s outer boundary does not substantially deviate 
from a circle. 

In chapter VI, the FE-BI formulation was developed for axially symmetric struc- 
tures. The finite element method for this problem was based on the CAP equations 
for generating the FE matrix system, and a boundary integral equation was used for 
the boundary condition on an arbitrarily shaped mesh termination boundary. Conse- 
quently the boundary integrals were no longer convolutional and a storage reduction 
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was, instead, achieved by exploiting certain symmetry properties of the boundary 
integral subsystem. Also, a new method was presented for circumventing the singu- 
larity problem associated with the CAP equations in the finite element portion of 
the system. 

A method was developed in chapter V for eliminating the internal resonance cor- 
ruption of the solution, a difficulty found with most boundary integral formulations. 
These resonances correspond to the eigenvalues of the integral operator, and may 
be also viewed as the cut-off frequencies of a resonator with conducting walls of the 
same shape as the mesh termination boundary. Because the eigenvalues become very 
closely spaced for electrically large structures, any scattering computations become 
unreliable unless the resonances are suppressed. 

In summary, specific contributions of this work include the development and 
implementation of the FE-BI for rectangular, circular and ogival boundaries as de- 
scribed in chapters III and IV in a manner taking advantage of storage reduction 
schemes. A new method was also presented for suppressing resonance corruptions 
existing in almost all implementations employing some form of a boundary integral 
equation over a closed surface or contour. This method involved the introduction of 
a complex wavenumber and was demonstrated for both 2-D and axially symmetric 
bodies. The implementation of the FE-BI method for axially symmetric structures 
as described in chapter VI was for the most part new and incorporated a scheme 
for treating the line singularity associated with the CAP equations. This resulted in 
real matrices for lossless scatterers, a property consistent with 2-D and 3-D imple- 
mentations of the finite element method. 
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7.2 Future Work 

The FE-BI technique presented in this work remains a promising technique for 
computing the scattering from ever larger structures. However, as the structure 
increases in size, so does the condition number of the resulting system due in part 
to its size. To reduce the number of unknowns, higher order basis functions must 
be used. For electromagnetic scattering, bases formed with hierarchical functions [5] 
show promise. Instead of standard functions, in which each undetermined coefficient 
corresponds to an unknown field quantity, those associated with the hierarchical 
function are comprised of both field quantities and field errors. For instance, a 
quadratic hierarchical function for a one-dimensional element takes the form 

^^JV‘(i)^ + a/(i) (7.1) 

i=i 

where f(x) is an arbitrary quadratic function, and the coefficient a is proportion to 
the deviation of the linear approximation to the quadratic one. Clearly, as higher 
order terms are added to the sum in (7.1), the smaller the coefficients become. Fur- 
thermore, the quadratic function f(x) may be chosen to be approximately orthogonal 
to polynomials of different order. Thus, in a finite element implementation, the off 
diagonal elements may become very small increasing the diagonal dominance of the 
matrix and leading to better matrix conditioning. 

The above idea can be extended to include functions other than polynomials. The 
physical optics solution could be employed as the fundamental solution and linear (or 
higher order polynomials) functions could then be employed as correcting functions. 
That is, the unknown coefficients in the solution would then be the deviation of the 
physical optics solution from the actual one. This could consequently aid in reducing 
the number of unknowns, particularly if the original approximation was reasonably 
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good. Note that this technique is applicable for the discretization of the boundary 
integral and finite element equations. 
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